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CONVERSION FACTORS, U.S. CUSTOMARY TO METRIC (SI) UNITS OF MEASUREMENT 


U.S. customary units of measurement used in this report can be converted to 
metric (SI) units as follows: 


EE 


Multiply 


by 


To obtain 


inches 25.4 - millimeters — 
2254 centimeters 
Square inches 6.452 Square centimeters 
cubic inches 16.39 cubic centimeters 
feet 30.48 centimeters 
0.3048 meters 
Square feet 0.0929 Square meters 
cubic feet 0.0283 cubic meters 
yards 0.9144 meters 
Square yards 0.836 square meters 
cubic yards 0.7646 cubic meters 
miles 1.6093 kilometers 
square miles 259.0 hectares 
knots 1.852 kilometers per hour 
acres 0.4047 hectares 
foot-pounds 1.3558 newton meters 


1.0197 x 1072 


kilograms per square centimeter 


millibars 
ounces 28.35 grams 
pounds 453.6 grams 
0.4536 kilograms 
ton, long 1.0160 metric tons 
ton, short 0.9072 metric tons 
degrees (angle) 0.01745 radians 
Fahrenheit degrees 5/9 Celsius degrees or Kelvins! 


lt obtain Celsius (C) temperature readings from Fahrenheit (F) readings, 
use formula: C = (5/9) (F -32). 
To obtain Kelvin (K) readings, use formula: K = (5/9) (F -32) + 273.15. 


A NUMERICAL MODEL TO SIMULATE SEDIMENT TRANSPORT 
IN THE VICINITY OF COASTAL STRUCTURES 


Marc Perlin and Robert G. Dean 
I. INTRODUCTION 
‘lee General. 


The need for reliable predictions of shoreline response to man-made or 
natural modifications is increasing due to environmental concerns and the 
rising cost of remediai measures. The capability of numerical modeling in 
addressing problems of shoreline response has advanced with improvements in 
wave climatology, programs to better understand sediment transport 
relationships, and improvements in numerical modeling. In-situ and remote 
sensing technology for the measurement of directional wave characteristics 
has been developed and applied, primarily within the last two decades. In 
addition to providing the necessary climatology, the resulting measurements 
have provided the basis for evaluation and refinement of directional wave 
prediction procedures. Studies such as the Channel Islands Harbor Longshore 
Sand Transport Study (Bruno, et al., 1981) and the Nearshore Sediment 
Transport Study (NSTS) (Gable, 1979) have yielded a better understanding of 
surf zone dynamics and the resulting sediment transport. The increased 
capacities of large computers and reduced computing costs combined with 
improved numerical modeling algorithms have resulted in an extremely 
promising potential for the numerical modeling of shoreline problems. 


Although it is doubtful that numerical modeling will ever replace 
completely the use of movable-bed physical models, the former type offers 
many advantages. The modeling of shoreline response is somewhat analogous to 
the problem of simulating storm surges in the coastal zone in which the scale 
effects and measurement difficulties essentially preclude physical modeling. 
For shorelines, the scale effects inherent in modeling sediment are well 
recognized and the costs of representing a substantial length of shoreline 
may be prohibitive. The laboratory representation of a realistic wave 
climate is at the forefront of technology. 


The investigation of shoreline response can best proceed by several 
approaches, with each approach selected for the particular strengths which it 
offers. Field programs are costly, usually because of the considerable 
equipment and the extensive time required, but these programs are essential 
for quantifying the values of constants or parameters, the forms of which may 
be available from laboratory measurements or theoretical considerations. 
Laboratory studies occupy a special niche by allowing the wave conditions and 
independent variables to be controlled readily, experiments to be repeated, 
and selected measurements to be conducted. Although, as noted before, scale 
effects are present in laboratory measurements of sediment transport, the 
physics governing the process should be the same. However, the relative 
magnitudes of suspended versus bedload transport in the laboratory and field 
may differ. Laboratory studies can also provide an excellent base for 
evaluating certain aspects of a numerical model, including wave refraction 
and diffraction and the resulting shoreline patterns due to, for example, the 
placement of a littoral barrier. Numerical modeling offers the capability to 


incorporate all the hydrodynamic wave-surf zone and sediment transport 
knowledge that is available from laboratory and field studies. Numerical 
modeling has the potential of providing accurate predictions of shoreline 
response to various structural and nourishment alternatives. Additionally, 
the possibility exists of employing numerical models and available field 
measurements to learn more about sediment transport mechanisms. In this 
latter mode, various candidate mechanisms or coefficients would be evaluated 
by determining the best match between measured and predicted shorelines and 
the bathymetry. Generally, this mode would require high-quality measurements 
of the forcing function (waves and nonwave-related currents) and the 
associated response (sediments) as well as the knowledge of appropriate 
conditions at the boundaries of the model. 


The present report documents the development and application of an 
n-line numerical model to investigate bathymetric response to time-varying 
wave conditions and shoreline modification. The model includes both 
longshore and onshore-offshore sediment transport. Based on laboratory 
results, a new distribution of longshore sediment transport across the surf 
zone is used. The wave climate is specified on the model boundaries which do 
not need to extend to deep water. Efficient algorithms are employed for 
representing wave refraction and diffraction. The equation of sediment 
continuity and transport are solved by a completely implicit algorithm which 
allows a large time-step. Specified sediment transport values or specified 
contour positions can be accommodated at the model boundaries. The model is 
suitable for investigating the shoreline response to a variety of 
modifications such as one or more groins, terminal structures, structures 
with variable permeability, and beach nourishment with or without terminal 
structures. 


2. Study Objectives. 

The objectives of the present study include (a) the documentation of 
state-of-the-art models, (b) the development and documentation of an improved 
model which includes the capability to represent n-contour lines and (c) the 
application of the model to several relevant coastal engineering problems. 


II. BACKGROUND 
This discussion describes significant contributions which either address 
numerical modeling of shorelines directly or provide improved capability for 
modeling. 
1. Wave Refraction (Noda, 1972). 
Noda developed an algorithm for solving the following steady state 
equation for wave refraction 


vxk=0 (1) 


in which ¥, the horizontal vector differential operator, and k, the wave 
number, are defined in terms of their components as 


Vos ie fg ee (2) 


k =ik +jk (3) 


where 7 and j are the unit vectors in the x and y directions respectively. 
Equation (1) can be expressed as 


a(ksine) _ a(kcose) 


in which e is the direction of the vector wave number relative to the x-axis 
and k denotes |k|. Noda expanded Equation (4) to the following form 


k cose = + sine ee = -k sine . + cose * (5) 
Since ok and oe are known from the angular frequency o, the water depth 
h, and the dispersion equation 
s@ = g k tanh kh (6) 


Equation (5) can be solved numerically, although there are problems of 
directional stability. The primary advantage of Equation (5) is that it 
allows the wave direction 6 to be determined on a specified grid, compared to 
unspecified locations that would be obtained by, for example, wave ray 
tracing. 


2. Crenulate Bays (LeBlond, 1972). 


LeBlond attempted to model the evaluation of an initially straight 
shoreline between two headlands into a crenulate bay. The model constitutes 
a one-line (shoreline) representation. The transport equation employed 
related the total sediment transport to total water transport in the surf 
zone as predicted by the formulation provided by Longuet-Higgins (1970). The. 
initial shoreline patterns resemble crenulate bays in nature; however, the 
predictions were found to be unstable for reasonably long periods of computa- 
tional time and did not approach a realistic planform. 


3. Crenulate Bays (Rea and Komar, 1975). 


Rea and Komar employed a rather ingenious system of orthogonal grid 
cells to provide a cell which locally is displaced perpendicular to the 
general shoreline orientation. A one-line representation was employed. A 
simple and approximate representation of wave diffraction was employed. 
Although the model yielded reasonable results for the examples presented, the 
unique coordinate system would not be suitable for a general model as the 
coordinate system must be "tailored" to some degree to conform to the 
expected shoreline configurations. 


4. General One-line Shoreline Model (Price, Tomlinson, and Willis, 1Y/c). 


Price, Tomlinson, and Willis' formulation consists of the sediment 
continuity equation and the total sediment transport equation 


i 0.70 a (nC) sina, COSa, 


Fie eT OLS ST = hoe (7) 


A) S 


in which E represents the wave energy density, (nC) the group velocity, 
a the angle between the breaking wave front and the shoreline, y, the 
Specific weight of water, p the in-place sediment porosity, and S. the 
specific gravity of the sediment relative to the water in which it is 
immersed. The subscript "b" represents values at breaking. 


Two formulations were presented by Price, Tomlimson, and Willis (1972). 
In the first, Equation (7) was substituted into the continuity equation and 
the results cast into a finite-difference form. In the second, the two 
equations were employed separately. The latter formulation was selected due 
to its simplicity and used for the results presented. 


Computations were carried out for the case of beach response due to the 
placement of a long impermeable barrier. The total sediment transport 
equation by Komar (1969) was used and the planform was calculated at 
successive times. Refraction was apparently not accounted for in the 
numerical model. To verify the computations, a physical model study was 
carried out for the same conditions using crushed coal as the modeling 
material. The comparison was interpreted as good for up to 3 hours; however, 
for greater times, substantial differences occurred and these were inter- 
preted as being due to wave refraction not being represented. The crushed 
coal was supplied to the model at the updrift end at a rate based on the 
Komar equation, and the results were interpreted as substantiating this 
relationship. However, the updrift end of the model beach receded substan- 
tially both in the numerical and physical models. In the physical model, 
this can only be interpreted as due to the Komar equation predicitions being 
less than the actual transport rate, possibly due to the low specific gravity 
(1.35) of the crushed coal. The predicted recession of the updrift beach is 
puzzling, although it could be due to a problem in properly representing the 
updrift boundary condition. 


Other one-line models for shoreline changes in the vicinity of coastal 
structures were developed by LeMehaute and Soldate (1977) and Perlin (1978). 
Perlin also developed a two-line model formulation, with one-line represent- 
ing the shoreline and the second the offshore. Dragos (1981) developed an 
n-line nodel for bathymetric changes due to the presence of a littoral 
barrier. 


III. THE NUMERICAL MODEL 


1. Description. 


There are several methods of modeling bathymetric changes due to the 
presence of a littoral barrier. An attempt can be made to either model the 
complete hydrodynamics and the resulting sediment transport or model using a 
combination of analytical and empirical sediment transport equations. The 
second method was chosen due to past relative success. 


At least two methods of employing sediment transport equations exist: a 
fixed longshore and cross-shore grid system where the depth is allowed to 
vary or a fixed longshore and depth system where the cross-shore distance is 
allowed to change. Although it may seem somewhat awkward, the latter system 
was chosen for the model. This method allows the modeler to think of 
bathymetric changes due to a littoral barrier in terms of the effect on the 
contours; i.e., the contour realinement due to the structure's presence is 
observed. One limitation of this approach, at least as it was applied here, 
is that each depth contour must be single-valued; it is not possible to 
represent bars. 


The next step in formulating the model was choosing the specific 
representation of the bathymetry. The model is an n-line representation of 
the surf zone in which the longshore direction x is divided into equal 
segments each ax in length. The bathymetry is represented by n-contour 
lines, each a specified depth, which change in offshore location according to 
the equation of continuity. There are two components of sediment transport 
at each of the contour lines, a longshore component, Qx, and an offshore 
component, Qy. Figure 1 is a definition sketch showing the beach profile 
representation in a series of steps and the planform profile representation 
and notations used. 


Implementation of the sediment transport equations requires knowledge of 
the wave field and the equilibrium offshore profile. A discussion of the 
refraction and diffraction schemes follows. The equilibrium profile is 
introduced when it is convenient. As an introduction to the logic used in 
the numerical model, a flow chart is presented in Figure 2. 


2. Refraction. 


A refraction scheme compatible with variable ay's was required because 
of the variable distance to fixed depth contours (as opposed to the more 
usual fixed grid system where a grid center has a longshore and offshore 
coordinate with a variable depth). One of the benefits of the n-line model 
is the ease with which the response of the contours to a particular wave and 
structure condition can be visualized. A fixed grid system and an 
interpolation scheme could have been used to obtain the wave field; however, 
this would have reduced accuracy and increased computation time. The scheme 
developed also saves computation time because it does not use differential 
products terms. 
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The first of the governing equations used is the conservation of waves 
equation 


+i, xk=0 (8) 


> > 

where vy is the, horizontal differential operator equal to i(a/ax) + jla/ay) 
in which i and j are the unit vectors in the x and y directions, respec- 
tively, and x is the longshore direction, with positive to the right when 
facing the water, y the offshore direction, with positive seaward, and z the 
vertical coordinate, with positive defined as upwards. For the steady-state 
case, equation (8) yields 


p) t) = 
# (ky) - & (k,) = 0 (9) 


where ky and ky are the wave number projections in the respective directions. 
Defining as the angle k makes with the y-axis positive in the counter- 
clockwise direction, the equation can be written in final form as 

3 3 : 

= (k cos @) = 5 (k sin @) (10) 
where @ = a + w (in radians). Noda (1972) and others have developed 
numerical solutions to expanded forms of equation (10). In the present 
study, equation (10) was initially central-differenced in the x-direction and 
forward-differenced in the y-direction with Snell's law used to specify the 
boundary conditions on the offshore boundary and one of the sides (i.e., the 
side of the wave angle approach). However, a numerical problem arose. The 
argument of the arcsine exceeded + 1.0 for large ay/ax. To overcome this 
problem, a dissipative interface was used on the forward-difference term 
(after Abbott, 1979). The final finite-differenced form of equation (10) is 
n+] ral 1 


= sin — ft sin e) 


8. 
led RT 


tele (1-2t)(k sin e). SL tan 


. A atl 
+ Tk Usain OV Fas - st ((k cose) 7 GF (k 7 had 


where t has been taken as 0.25. The past of j and the present of j wave 
angles are numerically averaged to give the 6},;. Newton's method is used 
to compute the wave number via the linear wave theory dispersion relation. 
In addition, numerical smoothing is used at the conclusion of the wave field 
calculation. This approximates in an ad hoc manner diffractive effects 
(lateral transfer of wave energy along the wave) which exist in nature but 
have been omitted due to use of the equation for refraction (equation 8). 
The smoothing routine is 


(12) 


The second governing equation used in the refraction scheme is 
conservation of energy. Neglecting dissipation of energy due to friction, 
percolation, and turbulence, this equation can be expressed as 


Welle ce) = 0 (13) 


where E is the average energy per unit surface area and Ce the group velocity 
of the wave train. Performing the operation indicated and replacing C by 
its components (Cosin e) and (C.cos e) results in the following: 


3 : 3 A 
a? (E Co sin e) + ay (EXGnZeose)e="0 (14) 


Assuming linear theory, 


E = agit (15a) 


where p is the mass density of water, g the, gravitational constant, and H 
the wave height. Dividing the equation by o> finite-differencing and 
weighting the forward-differenced term as before, and solving for the wave 
height, results in the following: 


r 


n+1 = i 2 2 
ig (C,cose). j (t)(H Cqcose) sy say + (1-27)(H Cqcose) 4) 


(15b) 


2 aY cine ¢3 “csi ae 


This equation is also solved by iterative techniques and the HEU and Hi ; 
are averaged at the conclusion of each iteration. oJ J 
Ce is determined by the linear wave theory relationship 
AG 2kh 
8 8 (ag sinh OKh) (16) 


where h is the water depth, k the wave number, and C the wave celerity. Wave 
height boundary conditions are input along the same boundaries as the wave 
angles using linear theory shoaling and refraction coefficients. The e's 
have been previously determined. In both equations (11) and (15) for a 
variable grid system, the points (i + 1, j) and (i - 1, j) need to be 
determined (i.e., because the y coordinates are not fixed, adjacent values 
with the same subscripts can be farther or closer to shore, therefore 
interpolation must be used). The actual values are found by searching the 

(i + 1) and (i - 1) cross-shore lines, finding the adjacent values in the 
positive and negative y-direction, and interpolating to determine the value. 


3. Diffraction. 


The diffraction solution (in the lee of the structure) used in the model 
is based on the method of Penny and Price (1952). Assumptions used in this 
method include a semi-infinite breakwater, which is infinitesimally thin, 
linear wave theory and constant depth. A definition sketch for wave 
diffraction is shown in Figure 3. The quantity THETAO represents the angle 
of wave incidence relative to the jetty axis, ANGLE represents the angle from 
the jetty at the point where the diffraction coefficient is to be computed, 
and RAD is the radial distance. The radial distance is then cast into a 
dimensionless parameter, RHOND (= 2n RAD/L), where L is the wavelength. This 
is equivalent to multiplying the radial distance by the wave number k. ~ 


The diffraction coefficient AMP is expressed as the modulus of the 
diffracted wave 


AMP = (Sum 1)2 + (Sum 2)@ (17) 
where 

Sum 1 = [cos (RHOND (cos (ANGLE-THETAO) )) (5 (1.0 + C. + S))] + 
[sin (RHOND (cos (ANGLE-THETAO))) . (- ; (S - C.))] + 
[cos (RHOND (cos (ANGLE+THETAO) )) peito ace siSiilce 
{sin (RHOND (cos (ANGLE+THETAO) )) G - (S-Ce))] (18) 

Sum 2 = [cos (RHOND (cos (ANGLE-THETAO))) . (- ; (S - C.))] + 
[sin (RHOND (cos (ANGLE-THETAO))) . (5 (1.0 +C, +S))]+ 
[cos (RHOND (cos (ANGLE+THETAO))) . (- : (S- C.))] + 
[sin (RHOND (cos (ANGLE+THETAO))) . (5 (1.0+Cp+S))] (19) 


In Equations (18) and (19), Cre and S represent Fresnel integrals and are 
computed in the model by means of an approximation after Abramowitz and 
Stegun (1965). 


Having obtained AMP, the wave height at the location in question is 
simply the product of the specified partially refracted incident wave height 
and AMP. The angle of the wave crest is computed assuming a circular wave 
front along any radial; this angle is then refracted using Snell's law. 


Throughout the refraction and diffraction schemes, the local wave 
heights are limited by the value, 0.78 x depth. 


Incident wave ray 


THETAO 


Location where diffraction coefficient, 
AMP, is to be calculated 


Ws 


Semi-infinite jetty 


Figure 3. Definition sketch for wave diffraction. 
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4, Sand Transport Model. 


a. Governing Equations. Three basic equations are used to simulate the 
sediment transport and bathymetry changes according to the wave field. The 
equation of continuity 


aq aq 
EN sie te Noy EL 
atu ay 0 (20) 


requires as input, knowledge of the longshore and cross-shore components of 
sediment transport. The total transport alongsnore has been measured by 
several investigators and many equations exist; however, the distribution of 
the transport across the surf zone is not well known. Fulford (1982) based 
on laboratory data from Savage (1959), developed a distribution of longshore 
sediment transport across the surf zone for the case of straight and parallel 
contours. Fulford's use of Savages experiment was based on two assumptions: 
1) the structure must be a total littoral barrier and 2) onshore-offshore 
sediment transport could be neglected. Test 5-57 was chosen because the two 
criteria were nearly met. Savage reported that the groin acted as a total 
littoral barrier for the first 35 hours of the test (i.e., no bypassing 
occurred prior to 35 hours). This does not mean that no onshore-offshore 
transport occurred because as the profile steepens on the updrift side, 
onshore-offshore transport does occur. However, it was assumed to be 
negligible. In addition, the initial profile had been molded to an 
equilibrium profile via 150 hours of waves. Thus, the two criteria required 
to develop an inferred longshore distribution of sediment transport were 
nearly satisfied. This distribution is shown as a dashline in Figure 4. The 
smaller "maximum" is believed to be an extraneous effect of a groin downdrift 
from the location in the experiment where the data were taken. Therefore, 
this feature was replaced by a monotonically decreasing, smooth curve as 
shown by the "altered" curve. To analytically represent this distribution, a 
function of the following form was chosen 


n 
a. (y= (8) Gam) en? (21) 


This type of equation is convenient because it is easily integrable, and by 
properly choosing the constant, B, the integral of the equation from zero to 
infinity can be required to equal a particular value. This too is highly 
desirable because, aS was done in the model, the integral is set equal to one 
and then multiplying by the value of the well-known longshore transport 
equation, the value of the transport at any location across the surf zone can 
be determined. Further investigation suggested a value of n = 3 to produce a 
curve similar to Fulford's curve. A more general form of the equation which 
allows more flexibility and curve fitting is 


Pe ia | Aaa : 
q(y) = Bly +a)" e CYp (22) 
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Figure 4. Distribution of sediment transport across the surf zone. 


where yp = distance to the point of breaking 


a = constant to allow sediment transport above mean water line (MWL) 
(swash transport or transport in region of wave setup) to be 
represented 

c = a constant establishing the width of the curve (to be determined) 


Ba = (causes fo (y) dy = 1.0) 
cy, 


O 


Based on Fulford's (1982) results and considering a to be proportional 
to the breaking height divided by the beach slope, the constant of propor- 
tionality was determined to be unity; i.e., a = hp/(ah/ay). Using equation 
(22) and a digitized version of the curve shown in Figure 4, a nonlinear least 
squares regression was carried out to determine the value of c. A Taylor's 
series expansion of the form 


k+1 (c,y) = £*(c,y) fesce AG (23) 


f 
where k and k + 1 represent the number of the iteration carried out. Least 
Squares regression minimizes the square of the difference between observed 
and predicted values with respect to a change in the parameter being 
computed, or 


2 
n= 


where fogs represents the observed values, which in this case is 
ax(y)ops- Carrying out the differentiation indicated and manipulating 
terms, Ac can be solved in terms of known quantities. 


An iterative procedure was then used by updating the values of 
fk(c,y), af/ac, and c until an acceptably small change inc results. For 
the data herein, the value of c was determined to be 1.25. The final form of 
sediment transport of a y location in the surf zone results for a shoreline 
with straight and parallel contours, as 


3 
(1.25)5(yp)° 
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This equation, which is also presented in Figure 4, predicts the relative 
transport at point y. To obtain the fraction of transport between two y 
coordinates, the integral of equation (25), from y; to y2, must be 

used. 


Ha sie Ly. = al/(t)25 9 ue 
Oe aa) = quiyildy =e) % Yiae ts 


o/ y 
1 1 3 


we ~El¥g# a)/(1.25 y,)] (26) 


Q,[ND] is dimensionless; therefore, to compute a value in, say, cubic feet per 
second, it must be multipled by the total transport along a perpendicular to 
the shoreline obtained from the total longshore transport equation used in 

the model 


5/2 ; 


0) Hp sin (2 ay) (27) 


See Appendix A for a discussion of the constant C'. It is noted that the 
transformation of qx(y) to qx(h) can be effected by multiplying by the 
one-dimensional Jacobian (ay/Ah). This latter form (qx(hn)) is more useful 
here because the present model simulates the changes in contour position (ay) 
rather than changes by depth (ah). 


In the numerical model, Qy (I,J) (see Fig. 1) is determined using 
equation (26) except for the shoreline contour, J=1, and the farthest 
offshore contour simulated, J = JMAX. The shoreline contour longshore 
transport, Qy (1,1), in order to include swash transport, uses equation 
(16); however, the first term is set equal to 1.0. The seawardmost contour 
transport, Qy (1,JMAX), in order to include any longshore transport not yet 
accounted for, neglects the second term of equation (26) (i.e., it accounts 
for transport from y(I,JMAX) to infinity). The dimensionless numbers are 
then multiplied by Q determined from equation (27). This method is based on 
parallel contours which may not exist. In order to compensate for the 
nonparallel nature of the contours (note that refraction does account for it 
as far as the wave field is concerned), the term sin (2ap) of equation (27) 
is replaced by sin (2a) shoreward of the breakpoint, where a represents 
the angle between the "local" wave angle and the "local" contour. It can be 
argued that for a spilling breaker, the remaining surf zone at any point 
"sees" a total transport similar to equation (27), where ap and Hp are 
the local values. The problem is that the constant of proportionality was 
determined for the entire surf zone and for nearly straight and parallel 
contours. This not being the case, the equation was altered on intuitive 
grounds to reflect the fact that the contours are no longer straight and 
parallel. 


2 | 


The second input required by the continuity equation to predict the 
bathymetric changes is the cross-shore sediment transport. The governing 
equation for onshore-offshore transport (after Bakker, 1968) is 


Q = ax C I St eTAV re ee PON (28) 
Vs OE i,j-l Tse] EQS 


where Co is an activity factor (inside,the surf zone = 10e2 feet per second 
for the OF Fototype simulation herein, 10 © feet per second for the physical 
model simulation) (see App. A. for a een pa We (i,j) is the 

positive equilibrium profile distance between y(i sneiy i,j-1), determined 
from the equilibrium profile used in the numerical pee ae ay2/3 (Dean, 
1977). See Appendix A for discussion of the value of A. The physical 
interpretation of equation (28) is that as this profile steepens (flattens), 
sediment is transported offshore (onshore). 


b. Methods of Solution. Three separate finite-difference techniques 
were used to solve the equations: 


(1) Explicit longshore-continuity and explicit cross-shore 
continuity; 


(2) Implicit longshore-continuity and explicit cross-shore 
continuity for half a time-step then vice versa; and 


(3) Implicit longshore-cross-shore continuity. 


An explicit formulation was first developed which used the refraction 
scheme, the distribution of longshore sediment transport across the surf 
zone, and the onshore-offshore sediment transport equation. Problems in 
addition to the usual ones which are encountered with explicit methods (e.g., 
computation time and cost) were immediately realized. In the explicit 
method, both transport computations are based on the former values of the 
contour locations and are completely uncoupled. Stability of an explicit 
scheme requires a small time-step. In addition, the noncoupled nature of the 
equations, in some cases, resulted in crossing of the contours due to the 
transport computed. 


It is logical to assume that an implicit formulation of the longshore 
transport equation used as input to the continuity equation along with the 
explicit onshore-offshore transport component would help the numerical 
stability (on the other half time-step, the longshore component would be 
computed explicitly and the onshore-offshore transport equation would be 
solved implicitly with the continuity equation). Although this scheme would 
be superior to the explicit procedure, it still would be susceptible to 
crossing contours. It should be noted that the magnitude of the coefficient 
used in the onshore-offshore equation is very important to the extent that 
the simulation models natural phenomena. If the coefficient is very small or 
vanishes, sediment will not move offshore and contours will cross because of 
the variation in the distribution of longshore sediment transport across the 
surf zone. If the coefficient is too large, the onshore-offshore transport, 
may become large enough that on a particular time step, an offshore contour 
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would move too far shoreward, thereby crossing an inshore contour or vice 
versa. Once the contours cross, not only does the bathymetry become 
unrealistic, but mathematically, the equation which computes the longshore 
distribution across the surf zone changes signs at some locations and the 
entire model becomes physically unrealistic. 


To circumvent these problems, an implicit scheme that simultaneously 
solves the three governing equations, was developed. Utilizing equation (26), 
and the one-dimensional Jacobian (ay/ah) to convert to Qy(h), the total 
longshore transport equation (27), the following equation is obtained, 


3/2 Sh203 Hie. 3/223 
(hy jen Hy A (ins Niges Hy, A 
Pane | ler (=~) a | 
Loa : b. ; b. 
x (5? | x sin (20 - 2a,) (29) 
1,J 


Qx(i,j) represents the sediment transport between depths h(i,j) and h(i,j-1) 
(see Fig. 1). The term in brackets represents the normalized distribution of 
longshore transport between h(i,j) and h(i,j-1); © is the averaged wave angle 
at the location of Qy(i,j) and ac is the local contour orientation angle. 
Defining everything except sin (20 - 2ac) as v(i,j) and using a superscript 
to denote a time step, this equation can be written 


gntl Save Suet ee 2aQ*t) 


Xe « Und) (30) 


The assumption has been made that the wave field (H and e) do not vary during 
the bathymetric changes over the time-step. Using the following trigonometric 
identities, 


Sime 2a eb)! = Sinyz2a.cos 2bs—-ncos. Za Sinseb (31a) 
cos 2a =H2\coselalel) (31b) 
sin 2a =i 2usinitanGoswa (31c) 


and recognizing that the following expression is an approximation 


1 nt+1 n+1 n n 
a iye al = hye s gaik Fee WA Bou (cea ee) 
Srna Negima ee sealed vee iey add Pela (32) 
3 (oy cs aaidlbin Ne Nae 
Wad) i-1,j 
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along with assuming that the change in the denominator is small for a 
reasonable time-step (the numerator has been averaged over the nth and n + 
ith time-steps), equation (30) results in 


n+1 n+1 nt n 
Oy j + (S3); 5 Vioj ra (S3); Yi-1,j os (RHS1) 5 (33) 
all 1 
where (S3). . = (5) (v. .) cos (2e) (2 cos ao.) —————————— 
Usd tod : 2 2\1/2 
(ox + AY 
n mR : 2 n n ; 
(RHS1); = (v; ;) (2 sin @ cos e)(cos a - ASS G AWG Aeshna al 


Here it has also been assumed that Cos ea does not change over the time 
step. Equation (33) is the final form of the longshore sediment transport 
equation prior to its use in conjunction with the other equations. 


Averaging y values on the nth and (n+1)th time-steps, equation (29) 
can be rewritten as 


oe 1 n+1 n n+1 n 
yg song a } (5. ie a Pies if i) 


where Const6(i,j) = Copfli,j). ax. This is the final form on the 


onshore-offshore sediment transport equation. 


The equation of continuity, finite-differenced for the nth and 
(n+1)th time-steps, can be written as 


ne ein 
Via) BO aE e 
ne [aati faa 1 nile sn n+1 n n+1 .,n n+1 n 
tated = ah gh art -0r Qh ag gna 
ne | MUS iad) Ateled Siete g ond ean dele 
(35) 


Defining Ry,j as 1/(2axah), inserting equations (33) and (34) into equation 
(35), and transferring all known quantities for the nth time-step to the 
right-hand side of the equation result in 


n+1 nt+1 nt+1 n+1 
Yi jt CAtRy 51535 Vag - COtRy 5)S35 4-1, 7 CPRE GS Sinn 5% 141, 
nt+1 il 


n+1 n+1 
+ (atR; 5 )S3 jeonst6; |) (; [ Ty ein ea i) 


1 er jon 


3 


1 nt+1 n+1 i 
ap (atR, ;Const6; 5.1) (; [ Vinj Fz Yi, j+l 1) = (AWARE); (36) 
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Equation (36) can be rewritten as 


(1+ U+ V+ Z1 + Z2) Tae - (u)yot? - (vjyntt 


i-1,J i+, 
— (ZFS y= (Zep sy = (AWARE), (37) 
where Ule= AtRs 59355 
PSOE Sia 
Z1 = (5) R; jConsté. | 
Z2 = (5*)R; jConst6, 41 


Equation (37) is a weighted, centered scheme in which yat is computed 
using a weighting of itself and its four adjacent grid "neighbors". The 
weighting factors (U, V, Zl, and Z2) are functions of the wave climate, the 
slope between contours, and the variables included in the original formula- 
tion. An investigation of a small gridded system demonstrated that by writing 
simultaneous equations, one for each yj,j, a banded matrix results. This 
matrix can be solved by LEQT1B, one of the available routines from the 
International Math and Statistics Library (IMSL). A schematic representation 
of the matrix A which results from the matrix equation [A][y] = [B] is 
presented in Figure 5. In this schematic, the large zeros represent 
triangular corner sections of all zeros and the 0...0 represents bands of 
zeros, the number of which is dependent on the number of contours simulated 
(the number of zero bands between either remote nonzero bands and the 
tridiagonal nonzero bands equals two less than the number of contours modeled 
(in both the upper and lower codiagonals of the matrix)). An inspection of 
the subscripts in equation (29) yields the reason the zero bands are required. 
The more j values (contours) used, the more y grids there are along any 
perpendicular to shore. This causes zeros to appear in the matrix between 
bands as the weighting factors await being used to operate on y"t!(i-1,3) 
and yn+l(i+1,j). For this reason, the expense of simulating an increasing 
number of contours is exponential. The LEQTIB routine, utilizes banded 
storage and saves both storage and computation time; however, the routine has 
no special way of handling the interior zero bands. One refinement which 
would save computation time would be to develop an algorithm to solve and 
store the matrix by taking advantage of these inner zero bands; however, it 
is beyond the scope of this project. 


Of course, the matrix requires boundary values on longshore extrem- 
ities and on both onshore and offshore boundaries. The longshore boundary 
conditions are treated by modeling a sufficient stretch of shoreline so that 
effects of a structure's presence are minimal. The y values along these 
boundaries can therefore be fixed at their initial locations. In the 
onshore-offshore direction, boundaries are treated quite differently. The 
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Note:. Size of matrix full storage mode 
[ (IMAX-2) (MAX) x (IMAX-2) (JMAX) ] 


Size of matrix panded storage mode 
[ (IMAX-2)(JMAX) x (20MAX + 1)]j 


Figure 5. Schematic representation of banded matrix if 
not stored in banded storage mode. 
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berm and beach face are assumed to move in conjunction with the shoreline 
position. The required sediment transport is then computed by the change in 
position of the shoreline. The two equations are 


Nitviaas He n 
Yi, 0 yi 0 a m Y¥; 11 (38a) 
Mtl ee Berm Ax n+] 


The offshore boundary is treated by keeping y"+](i,jmax) (the contour 
beyond the last simulated contour) fixed, until the angle of repose is exceeded. 
Then, the y"*t!(i,jmaxt+1) is reset (at the conclusion of the n + 1 time-step) 
to : position such that the slope equals the angle of repose. Note that 
yntl (1,0) is represented in the program by YZERO;. 


There are also no-flow boundary conditions required at each of the 
structures being modeled. These are imposed on the adjacent y-grid points which 
are located downdrift (i.e., in the shadow zone) of the structure and shoreward 
of the structures' seaward extremities. They are imposed by setting S3; j of 
equation (33) and DISTR; ,j (the term in square brackets in equation (29) 
equal to zero, thereby causing Qx(i,j) to be zero (i.e., the no-sediment flow 
condition). This boundary condition is imposed automatically for every 
shore-perpendicular structure. 


It was found that even with the implicit formulation, high frequency 
oscillations occurred in the y values immediately updrift and downdrift of the 
structure. The solution did not "blow up"; however, on larger time-steps 
"sloshing" (oscillating) did occur. Part of this problem was due to the 
boundary condition at the structure which had been such that either no sand was 
allowed along a contour line or the sand determined by the equations was allowed 
to be transported. Because of the very large angle which existed around the tip 
of the structure when a contour first exceeded the length of the structure, very 
large amounts of sediment transport were predicted. In the nature where analog 
sand transport rather than digitized transport occurs, this does not happen. 
Therefore, the boundary condition was altered to constantly allow sand transport 
around the end of the structure in proportion to that part of the contour 
representation which exceeded the structure (i.e., the transport was calculated 
for the location at tip of the structure as if the structure was not there and 
then a proportion of this value was allowed to bypass). Although the transport 
around the tip of the structure is based on the values from the past time-step, 
it more closely simulated the natural phenomenon. 


Additionally, a dissipative interface is used on the y values as follows: 


Yev= 2) yeas Ces alae siga (t) y (39) 


ieee) itl ,j 


where Tt was again taken as 0.25. It is noted that only high frequency 


oscillations in y are affected by the use of equation (39); the total sum 
of y values is not affected. Also, in all the dissipative interface 


(ATE 


schemes used, if a boundary point is being computed, either a forward- 
difference or a backward-difference of equation (39) is used (after 
Abbott, 1979): 


Backward: oy. . = (t)y,_) ay (1 - t)y (40a) 


lad i,j 


Forward: y =} (Ge Ne + (1 -t)y (40b) 


i,j itd isJ 


IV. SIMULATIONS AND VERIFICATION 


Several simulations were run; two were attempts at verifying the 
numerical model, the others were run to gain insight. Because a complete 
data set does not exist, only the available data are compared. The first 
modeling effort was to simulate the physical model tests of Savage (1959). A 
second set of cases was run for shore-perpendicular structures. Next, an 
effort was made to model sediment transport in the vicinity of a hypothetical 
dredge disposal site in the 1]1- to 14-foot depths off Oregon Inlet. Finally, 
the Channel Islands Harbor Longshore Transport Study (Bruno, et al., 1981) 
was modeled. Bathymetric changes were closely monitored during this study; 
however, the wave climate (H, e, T) used was determined from the Littoral 
Environmental Observation (LEO) data and uncertainties exist as to the 
accuracy of the data. 


1. Simulation of Savage's Physical Model Tests. 


The numerical model was used to simulate one of the physical model tests 
of Savage (1959). Test 5-57 was simulated numerically for a 10-hour period. 
In this physical model, the mean sediment size was 0.22 millimeters, the wave 
height averaged 0.25 feet, the wave period was 1.5 second, the wave angle was 
30° (at a depth of 2.3 feet), and the groin was approximately 9.5 feet from 
still water to its seaward limit. Corr was held constant at 10-4 feet 
per second throughout the profile for this simulation. The offshore profile 
is presented in Savage (1959). Figure 6 represents three of the eight 
contours simulated. Note that the initial 0.3- and 0.5-foot-depth contours, 
in the numerical representation are too far seaward by approximately 2 feet. 
This is due to the h = Ay2/3 equation as compared to the equilibrium 
physical model profile. Realizing this, it is the shape of the contour which 
must be used as an indication of the numerical model predictions. The 
general trend of the contours is similar, although the numerical model 
contours are displaced farther seaward as expected. The major differences 
are in the diffraction zone. 


2. Several Runs Using Shore Perpendicular Structures to Demonstrate Effects 
of Altering Some of the Pertinent Parameters. 


In the following simulations, the models were run until their 
near-equilibrium values were achieved. Coefficient Corr was not a function 
of depth (beyond the surf zone) but was held constant throughout the 
simulated area. Important variables are as shown in the figures. Only one 
wave condition (Ho = 3 feet, T = 7 seconds, and a deepwater wave angle ag 
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of 60°) was used as input for all four cases. Case 4.2a used an equilibrium 
shape factor A of 0.0899 and one groin. Case 4.2b was similar to 4.2a with 
the only modification being, that the A value was changed to 0.1486. In this 
way, a direct comparison was made based only on the shape of the equilibrium 
profile. Cases 4.2c and 4.2d used A-values of 0.0899 and 0.1486, 
respectively, but this time three shore-perpendicular, evenly spaced 
structures were simulated. 


a. Comparison of Cases 4.2a and 4.2b. The most obvious difference between 
Figures 7 and 8 is the volume of sand impounded updrift and eroded downdrift. 
This is due to blockage of more of the active transport zone in the second 
case (i.e€., a shorter groin is required for an equivalent performance on a 
steeper beach). The next obvious difference is the size of the perturbation 
which exists in the offshore contours. Clearly, case 4.2b is more perturbed 
and this is expected because larger offshore transports occur due to the 
steepening on the updrift side. Conversely, this means less sediment is 
initially bypassed (and along with the downdrift requirement for larger 
volumes of sand) causes larger erosional features in case 4.2b. Another 
interesting feature is the downdrift fillet which occurs in the third, 
fourth, and fifth contours. The fillet is due to the shape of the sixth 
contour which occurs because of the inability of the wave to transport more 
sediment (due to the reduction in wave height and angle in the diffraction 
shadow zone). The remaining difference is also due to the volume of sediment 
being impounded; i.e., the distance and extent of change the presence of the 
groin causes upcoast and downcoast. 


b. Comparison of Cases 4.2c and 4.2d. The variations between cases 4.2c and 
4.2d are very similar to the differences between cases 4.2a and 4.2b as would 
be expected with a groin field (here, three groins) as compared with a single 
groin (see Figs. 9 and 10). There is, however, one additional feature which 
can be attributed to the additional groins. Note that in the direction of 
littoral drift, the size of the fillet is decreasing. This is due to the 
updrift beach having an uninterrupted supply of sediment while the downdrift 
groin compartments are supplied sand at a rate determined by the bypassing. 
Part of this feature may also be due to the system not having attained 
complete equilibrium. 


The effects of the fixed boundary conditions are evident on all cases 
run. In these example cases, the boundaries are clearly too close to the 
structure to provide a proper representation of the fillet contours. 


Sic Simulations of Sediment Transport of Dredge Disposal in the Vicinity of 
Oregon Inlet. 


Hypothetical dredge disposal movement in the nearshore but beyond what 
is normally the surf zone at Oregon Inlet's adjacent beach to the south was 
modeled. In order to do these simulations, the program was altered such that 
for every nth iteration (time periods), the contours were shifted seaward 
to simulate the addition of dredged sediment disposal. The program presented 
in Appendix B does require slight modification to simulate this situation. 


In general, the fifth and sixth contours were shifted seaward on a 
monthly basis to simulate the disposal of 121,000 cubic yards of sediment. 
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In all these simulations, the following variables were held constant: (a) a 
time-step of 3 hours, (b) a shoreline length of 10,000 feet, (c) a longshore 
space-step of 200 feet, (d) an A value of 0.15 foot!/3 for the equilibrium 
profile (see Fig. 11), (e) a berm height of 5.3 feet with a beach face slope 
of 0.05, and (f) a duration of 1 year. The wave climate was provided by 

the U.S. Army Engineer Waterways Experiment Station Wave Information Study 
(WIS) 1975 data and was initiated at different times of the year as indi- 
cated in the specific cases below. All simulations, prior to any addition 
of sediment, used the bathymetry shown in Figure 12. The shoreline 

(relative to mean low water, MLW) was scaled from a bathymetry-topography 
survey provided by the U. S. Army Engineer District, Wilmington. The initial 
offshore bathymetry was computed according to the equilibrium profile and the 
0-foot contour; i.e., the profile was shifted seaward or landward, 
accordingly, (see App. C.) The boundary profiles were fixed throughout the 
simulations. The variation of COFF outside the surf zone was used because of 
the importance of the time rate of change in this simulation. Table l 
presents the percentage of sediment which moves out of the control volume 
(i.e., imaginary boundaries around the area where sediment was added) 
directly onshore and the percentage of sediment remaining in the control 
volume at the conclusion of the simulation for each of the cases. In 
addition, a seventh (case 3) and eighth (case 4) were modeled. In Case 3, 
the only difference was that sediment was placed at the 11- and 14-foot 
contours. Case 4, however, was quite different and will be described in 
detail later. It has a 20,000-foot shoreline, a longshore space-step of 400 
feet, and sediment was added on a weekly basis. Also, the resolution in the 
profile was better. 


diem opecitic Cases. 


(1) Case 2.a. In order to provide insight for the interpretation of the 
other modeling efforts, a simulation of the shoreline evolution using the 
January to December WIS time series, with no addition of sediment, was 
carried out. As expected, the contours almost attain an equilibrium planform 
shape (i.e., straight and parallel between the fixed end profiles; they do 
not, however, become aligned parallel to the base line because of the end 
conditions). Because of the scales involved, alongshore versus 
onshore-offshore, plotting the contours without distortion does not yield 
much information. Appendix C provides a listing of the final contours for 
all the cases modeled. 


(2) Case 2.b. The only difference between cases 2.a and 2.b is the 
suppression of the WIS wave angle which was set equal to zero (i.e., wave 
crest approach is shore-parallel at the offshore boundary of the model). 

This does not cause the longshore sediment transport to vanish completely. 
There are still local gradients in the contours which cause refraction and 
relative angles between wave crest and contour, thereby driving the longshore 
sediment transport (even if refraction was not considered, the local angle 
between the wave crest and contour would cause sediment transport). Note the 
larger onshore transport (Table 1) for this case compared with Case 2.a. 

This is due to the reduction in longshore transport caused by the wave angle 
of 0°. The model still tries to smooth the contour lines; however, more of 
the smoothing for the present case must be done by onshore-offshore transport. 
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Table 1. 


Summary of results at Oregon Inlet. 


Case Pct Onshore out of Pct Remaining 
No. Description control volume in control volume 


* After 17 weeks, the addition of sand caused contours to cross. 


WIS waves, Jan. 


No sediment added, 
WIS waves 
Jan. - Dec. 


No sediment added, 
WIS waves {a =0-) 
Jan. - Dec. 


121,000 yd3 added 
monthly, WIS waves 
Jan - Dec. 


121,000 yd3 added 
monthly, WIS waves 
Apr. - Mar. 


121,000 yd3 added 
monthly, WIS waves 
July - June. 


121,000 yd3 added 
monthly, WIS waves 
Oct. - Sept. 


121,000 yd3 added 
monthly at the 11- 
and 14-foot contours 


27,923 yd3 added 
weekly on the 7- 
8-, 9-, and 10-foot 
contours, WIS waves 

Jan. - Dec. 


sediment added was 363,000 yds. 


not rerun. 


- Dec. 


Onshore Movement 
(992 yd3) 


Onshore Movement 
(1624 yd3) 


Bila: 
(460,264 yd3) 


32.1 
(466,160 yd3) 


28.6 
(415,784 yd3) 


OD 
(395,556 yd3) 


8.9 * 
(32,164 yd3) 


19.0 
(275,796 yd3) 


Increase 
(14,148 yd3) 


Increase 
(9,356 yd3) 


38.6 
(559,984 yd3) 


36.9 
(535,392 yd3) 


47.0 
(682,088 yd3) 


46.8 
(670,848 yd3) 


78.0 
(283,016 yd3) 


47.4 
(687,525 yd3) 


Prior 


Problem was rectified; however, case was 
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(3) Case 2.cl. In this simulation, sediment is added to the system each 
month. It was simulated by advancing the 7- and 11-foot contours on a 
monthly basis to represent 121,000 cubic yards per month. Specifically, the 
sand volumes were "tapered" starting at the center of the nourished area over 
a distance of + 2,700 feet from the center. Table 2 presents the monthly ay 
values for the blocks between the 7- to 11-foot contours and the 11- to-14 
foot contours. Figure 13 shows the planform ay values added monthly. WIS 
waves were used with the sequence being the normal calendar year, January 
through December. 


eames 

ep = 4.0 ft 5 

increment 200 ft= Width of 
grid 


OFFSHORE 
DIRECTION 
145.8 ft 


I = 14 16 18 50 22 ZZ as mr 32 34 36 38 


Figure 13. Monthly incremental values of Ay due to dredge disposal 
illustrated for the block between 7- and 11-foot contours. 


The initial and final fifth and sixth contours have been plotted in 
Figures 14 and 15. The first figure has no distortion; the second is 
distorted 10 to 1. The simulation predicts that 31.7 percent of the dredge 
disposal will move shoreward out of the control volume. An additional 29.7 
percent efflux occurs in the offshore and longshore directions, leaving only 
38.6 percent of the total amount of sediment added remaining in the control 
volume. It is not clear what quantity of the sediment leaving in the 
longshore direction would reach shore. It is conceivable that most of this 
sediment would eventually reach the surf zone. The rate at which this 
material would move ashore would be expected to be slower than the rate at 
which the large mounds would move ashore because the deviation of the profile 
from equilibrium is much less. 


(4) Cases 2.c2, 2.c3, and 2.c4. The next three simulations were the 
same as 2.cl except the time series of wave events has been seasonally 
altered. Cases 2.c2, 2.c3 and 2.c4 use the 1975 wave climate from April 
through March, July through June, and October through September, respectively. 
The maximum variation is about 5 percent for the sediment volume moving 
onshore, and about 10 percent for the volume remaining. The variation in the 
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Table 2. Monthly values of Ay for the steps located between the 7- to 
10-foot contours and the 11- to 14-foot contours. 


Value of I Monthly ay value (ft) for steps between 


7- and 11-foot contours 11- and 14-foot contours 


26 
Zd5Cl/ 
24,28 
23,29 
22,30 
21531 
Z20R32 
19,33 
18, 34 
17535 
16,36 
152537 
14, 38 
13,39 

All Others 
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quantity moving onshore could be caused by waves that first tend to move more 
sediment longshore; then, the waves that transport more sediment onshore have 
a less out-of-equilibrium profile to cause movement upon. The variation in 
percentage remaining is due to the variation of the time series of the wave 
climate, with the last month in the simulation being especially important. 


(5) Case 3. Instead of extending the 7- and 11-foot contours monthly to 
simulate the disposal of dredged sediments, the 11- and 14-foot contours were 
extended (194.4 feet each at the center of the disposal area). This case was 
modeled because the larger available dredge could not dump in more shallow 
water. The reduction and increase in the percent of onshore volume and the 
percent volume remaining (8.9 percent and 78.0 percent, respectively) 
demonstrate the sensitivity of the depths investigated. Qualitatively, these 
depths are the depths to which offshore bars occur along the Atlantic U.S. 
coast. 


(6) Case 4. Further investigation of the disposal process demonstrated 
the need for an 11,000-foot shore-parallel disposal length with the sediment 
placed at the 11-foot contour building to about 7 feet. It was decided to 
model this physical situation also. The total shoreline length was changed 
to 20,000 feet, and the space step to 400 feet; the length of the disposal 
area in the longshore direction was increased to 10,800 feet. The resolution 
in the vicinity of the depths of the dump was improved by adding the 
additional contours and the profile is shown in Figure 16. As in the other 
seven cases, 1,452,000 cubic yards was added annually to the system; however, 
the addition was accomplished on a weekly basis (27,923 cubic yards per 
week). Sediment was still added by extending the contours seaward, but 
rather than placing one-fourth of the sediment at each of the four contours, 
the volumes were determined based on the trapezoidal cross section shown in 
Figure 17. This cross section more closely resembles the disposal mound 
formed by hopper dredging. The incremental values Figure 18 show, in 
planform, the extension of the contours to simulate the weekly sediment 
addition at the 8-foot contour. 


A schematic illustration of the sediment transported from the nourished 
region is presented in Appendix C. Nineteen percent of the sediment added 
moved directly onshore out of the control volume. 


b. Conclusions for the Movement of Disposed Sediment in the Vicinity of 
Oregon Inlet. The computer simulations, tempered with engineering judgment, 


demonstrate that between 15 and 35 percent of the material added to the 7- 
and 11-foot contours,or to the 7- 8- 9-, and 10-foot contours would be 
transported into the nearshore transport system during the first year. If 
the disposal process was continued, the system would approach steady state in 
terms of the volume of deposited material residing offshore. 


For the case of sediment addition at the 11- and 14-foot contours, the 
computer simulations, tempered with engineering judgment, show that between 5 
and 25 percent of the material added would be transported into the nearshore 
transport system during the first year. 
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28.05 pct. 9.0 ft 


10.0 ft 


VOL:35.16 pct. 11.0 ft Elevation 


Figure 17. Shore-perpendicular cross section of disposal 
mound. The volumes represent the volume per- 
centage of the trapezoidal section between contours 
and therefore, the quantity of sediment added to 
the 7-, 8-, 9-, and 10- foot contours. 


Contour 
Depth =| ft 
Increment 400 ft = width of grid OFFSHORE 


DIRECTION 
he 3 


Figure 18. Incremental values of Ay due to dredge disposal, 
illustrated for the block between 8- and 9-foot 
contours (case 4). 
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4. Simulation of the Longshore Sand Transport Study at Channel Islands Harbor, 


California. 


The CIH Longshore Sand Transport Study (Bruno, et al., 1981) was the 
only field study found suitable for verification purposes. Wave data 
collected included the LEO data and a two pressure-sensor gage array. 
Although the pressure gages were not in operation throughout the study, it 
was expected that the data they produced would be superior to that of the LEO 
data. However, these data were not available in a reduced form, so the LEO 
data were used. An adjustment of 11° was made to the breaker angle to orient 
. the angle with respect to the base line, rather than to the local shoreline 
orientation angle. Observations had been taken twice daily at three 
locations; the middle location was used (observer No. 5714). Waves which 
approached the shoreline at angles too large to have originated in a depth of 
10 meters, according to Snell's law, were set equal to 90° at that depth 
(crest of wave perpendicular to the baseline). The 10-meter depth was chosen 
because it is the approximate depth at the tip of the offshore breakwater 
(for this reason, it was also chosen as the depth of the step beyond the y(I, 
JMAX + 2)th contour). It was assumed that each of the two daily 
observations occurred for 12 hours and using a time-step of 6 hours, this 
meant two time-steps per wave. In cases where parts of the wave data (Hp, 
ab, or T) were missed by the observer or were equal to zero, the data were 
ignored (no computations were made), but the time was included. Because the 
time rate of change is important for this simulation, the variation of Corr 
outside the break point was used. 


The period chosen to model was 20 April through 1 December 1976. The 
initial survey was taken after dredging of the sediment trap and for this 
reason was known to be out of equilibrium. The bathymetric surveys were con- 
ducted using several methods, the most advanced being a Lighter Amphibious 
Resupply Cargo vessel (LARC) proceeding along shore-perpendicular lines 
(approximately in the vicinity of each survey station) taking fathometer 
readings every 10 seconds, with positioning systems trilaterating the 
vessel's position concurrently. These data were recorded on tape. The 
beach-face data were taken using standard surveying methods. Because the 
data fluctuated randomly about the stations, depending on the speed of 
the craft, the (x, y) coordinate positions had to be altered to fixed changes 
in x and y. This was accomplished using an interpolation routine. The x 
values were made to coincide with the stations used in the surveys, and the y 
values were determined at 100-foot intervals beginning from the base line. 
Stations 100+00 and 118+00 were located at the north jetty and termination of 
the detached breakwater, respectively (these correspond to I values of 16.5 
and 34.5 in the model). See Figure 19. 


Monotonic profiles of the form h = Aly - yde1)2/3 were fit to the data 
along each station line. "ydel" represents the zero location of the fitted 
shoreline, the value of which was unknown. Because dredging had been done in 
the lee of the breakwater, there was no reason to expect the A value to 
correspond to the value upcoast where the influence of the structure and the 
dredging was negligible. For this reason, the profiles of Stations 122+00 
through 134+00 were evaluated separately to determine an A value for the 
equilibrium profile to be used in the numerical model. For the detailed 
method used (LaGrange Multipliers and Newton-Raphson Method for nonlinear 
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Figure 19. Idealized numerical model representation of offshore 
breakwater at Channel Islands Harbor, California. 
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equations) and the computer programs see Appendix D. The two values obtained 
for the surveys of 20 April and 1 December 1976 were averaged to obtain the 
value used in the model, A = 0.2606. Stations 101+00 through 121+00 were 
treated separately for the purpose of obtaining values with which to initial- 
ize those parts of the contours in the model and for comparison of the model 
predictions with the prototype values. Note that although the breakwater 
extends only to about Station 118+00, the influence of the structure and 
dredging extends beyond that location and so, although arbitrary, the 121+00 
station was chosen as the dividing line. The initial and final values of the 
scaling parameter A for the profiles were 0.3233 and 0.3528, respectively. 
Because the initial shoreline is so irregular, a discontinuity between 121+00 
and 122+00 is not evident. 


One further idealization was made. The jetty-breakwater system was 
idealized as shown in Figure 19. This was required to simplify the physical 
situation, and although waves, currents, and sediment do pass through the 
opening in the prototype, it is hoped that they are of secondary importance. 


The results of the numerical modeling of Channel Islands Harbor are 
presented in Figures 20 and 21. The first figure presents the shoreline 
contour (depth = 0); the second figure presents the farthest offshore, 
modeled contour. In both cases, the initial shoreline represents the model 
and prototype (after fitting of the profiles). The initial shoreline contour 
is further offshore along the section of beach beyond the end of the 
breakwater, while in the lee of the breakwater, as would be expected after 
dredging, the shoreline is closer to the base line. The final prototype 
contour has undergone erosion along the reach beyond the tip of the 
structure, and accretion in the ‘ee. 


The model's shoreline contour has undergone similar changes, and on the 
average, represents the final prototype contour quite well. The JMAXth 
contour has been displaced quite similarly to the shoreline contour with 
shoreward movement (erosion) along the reach beyond the tip of the breakwater 
and seaward movement (accretion) within. It appears that the final model's 
shoreline has predicted too much erosion and not enough accretion. Several 
parameters could be incorrect, with the onshore-offshore sediment transport 
rate coefficient, Corff, perhaps the most likely. Overall, the model seemed 
to predict reasonable values of the contours. 


V. SUMMARY AND RECOMMENDATIONS 


Some of the parameters that the model does not include are important and 
should be mentioned. As stated previously, the model does not include bar 
formation. This is precluded by an n-line system. There are no provisions 
for water level fluctuations or currents. Improvement to the model could 
also be facilitated with better longshore and cross-shore sediment transport 
relationships. A more reliable equation for distribution of sediment 
transport across the surf zone would also be helpful (or further testing and 
calibration of the equation proposed herein). Finally, combining refraction 
and diffraction using equations to predict their combined effect would 
improve the wave field. The program was constructed such that improvement 
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Figure 21. CIH simulation of (JMAX)t" contour, 20 April - 
1 December 1976 (from LEO data). 
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could be accomplished with minimum effort. Therefore, if a more suitable 
equation becomes available, the change of a subroutine should be sufficient 
for implementation of the equation. 


Although the model is limited by the omission of the aforementioned 
parameters, it is reasonably correct. The ability to simulate various 
physical situations (shore-perpendicular structures, beach fills, breakwater 
and shore-perpendicular structures) has been demonstrated. In the CIH 
simulation where the data were first transformed to monotonically decreasing 
contours and LEO wave data were used, the model still predicts the prototype 
shoreline changes in a reasonable fashion. 


Further research and model development should include exercising the 
model in a number of different situations. Several theoretical cases should 
be simulated, which if analyzed properly, would provide a tool for the 
coastal engineer. Combined refraction and diffraction should be included, if 
possible, along with any of the aforementioned parameters which have been 
omitted and for which relationships exist. Perhaps the most difficult prob- 
lem to researchers working on modeling sediment transport in the vicinity of 
structures is the availability of field data. High-quality concurrent wave 
and bathymetric change data in the vicinity of coastal structures do not 
exist. One suggested field experiment is to monitor changes both updrift and 
downdrift of a jettied inlet which has a bypassing plant. Monitoring should 
begin immediately after bypassing, when the profiles are out of equilibrium. 
The recorded bathymetric and wave data would then provide data with which to 
calibrate, verify, and evaluate the existing models. 
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APPENDIX A 


DISCUSSION OF CONSTANTS AND SOME OF THE VARIABLES 
REQUIRED BY THE MODEL 


Establishing the grid-contour system requires several variables. IMAX 
represents the number of cross-shore grid lines desired and JMAX the number 
of contours simulated. DX represents the spacing between the IMAX grid lines 
and DY the spacing between the contours. DX is a value which must be chosen 
along with IMAX and JMAX such that sufficient detail is obtained where 
necessary (e.g., in the shadow zone, if diffraction effects are believed to 
be very important, DX must be assigned a sufficiently small value so that at 
least some points lie within the shadow zone for the larger wave angles). DY 
is not a constant, but a dimensional array which is computed by the model 
according to the contour location. Once the depths of contours to be modeled 
are chosen, the initialization of DY and the y values are computed with the 
following equation after Dean, 1977 


h=A y2/3 GAS) 


where h is the depth, y is the offshore distance and A is the scaling parameter 
Dean gives values for A for several diameter sediments; however, if long-term 
beach profiles are available for the site being modeled, the modeler may want 
to choose a slightly different A value to more closely match the site-specific 
beach profile. Figure A-1 presents values of A versus diameter (after Moore, 
1982). The model is programmed to input the h({I,J) values (depths as shown in 
Figure 1, called DEEP (I,J) in the program) read in the value of A (called 
ADEAN in program) and it then computes the y values. Also shown in Figure 1 
is the height of the berm (BERM) and this value, along with the beach-face 
slope (SFACE), is required as program input and can be obtained from beach 
profile site data. Because the model does not include water level 
fluctuations such as tides, all values are to be referenced to a chosen 

datum. Other geometrical constants depending on the site include SJETTY (the 
length of the jetty), MMAX (the number of structures to be input), and IJET 
(M), M = 1,2,...MMAX (the smaller I value adjacent to the mth structure's 
location). If no structure is required, as in a beach fill, the value of 
SJETTY must be entered as 0.0, with MMAX and IJET (M) entered as 1 and 
(IMAX/2), respectively. As set up presently, the groin locations must be 
equally spaced. 


One constant used throughout the program is the breaking wave criteria 
(CAPPA in the program) equal to 0.78. It is required in several different 
computations and always governs the maximum wave height allowed according to 
the depth. 


Another group of variables assigned values within the program is the 
sediment and fluid properties. These include fluid mass density, sediment 
mass density, porosity, and the angle of repose (e.g., RHO = 1.99, RHOS = 
5.14, POROS = 0.40, and REPOSE = 32°, respectively). The values can easily be 
changed to reflect site conditions. 
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Another very important set of constants is the constant chosen for the 
longshore and cross-shore components of sediment transport. Equation (27), 
the total longshore transport equation, contains the constant C' equal to 


ant yl/2 


Gu oy (A-2) 
ib. 2a) Se) Ge 


2 K 
Ps p 


where K = 0.77 (Komar and Inman, 1970) 
g is the acceleration of gravity (32.17 ft/sec2) 


pg and p are the mass densities of the sediment 
and the seawater (5.14 and 1.99 slugs per cubic feet, 
respectively 


p is the porosity (0.40), and 
K 7s taken as 0.78. 


Using these values to compute C' (TKSI in the program), a value of 0.325 is 
obtained. It is stressed that if any of these values are different for the 
site to be modeled, they should be changed and the program will compute 
another value for C’'. 


The parameter Corf is an “activity factor" which, based on earlier work 
primarily within the surf zone, was found to be 


eyed 
Corr = 10s sht/se h<h 


b 

To generalize this concept for transport seaward of the surf zone, the 

wave energy dissipation per unit volume was utilized as a measure of 
mobilization of the bottom sediment. Inside the surf zone, the dominant wave 
energy dissipation is caused by wave breaking; outside the surf zone, the 
dominant mode of wave energy dissipation is due to bottom friction. These 
two components will be denoted by Dy and Do, respectively. 


(a) Energy Dissipation by Wave Breaking. The wave energy dissipation 
per unit volume by wave breaking, Dj, is 
De epee (EsCR) (A-3) 
1 7h ay G 


which, employing the spilling breaker assumption (H = kh) within the surf 
zone, can be shown to be 
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5 g3/2 2yl/2 ah 


Bie fe a (A-4) 
or 

D, = < ogo! 2 Za3/2 (A-5) 
in which A is the scale parameter in the equilibrium beach profile 

h (y) = ay2/3 (A-6) 


(b) Energy Dissipation by Bottom Friction. The wave energy dissipation 
per unit volume due to bottom friction, Dy, is 


Asi Malt tes 3 
Do hy = ROC MY UE (A-7) 


in which C¢ is a bottom friction coefficient, up is the bottom water 
particle velocity and the overbar indicates a time average. For linear 
waves, equation (A-7) can be reduced to 


Deeneey One te (A-8) 
2 br h “f sinh kh 


The activity coefficient Corr: outside the surf zone, is expressed as 


1 2 uc 
Corr oa T 0. x 10 ft/s, h > hy (A-9) 
1 
Cao° 3 
paageieneiel f H x 10° (A-10) 
OFF 5r 3/2 2,3/2, \sinh kh 
GivtaicsAnaan 


in which fT is a parameter relating the efficiency with which breaking wave 
energy (which occurs primarily near the water surface) mobilizes the sediment 
bottom (0 < r <1). Herein, © is taken to be one. 


Figure A-2 presents an example of the variation of the activity 
coefficient versus relative depth for a particular wave period and deep water 
wave height. It is seen that the activity coefficient reduces rapidly with 
increasing depth. 


The value of Corr for the physical modeling of Savage's (1959) data 
was taken as 10-4 feet per second. Perlin (1978) presents some rationale 
for choosing a value of Corr; however, very little testing has been done 
and none is based on actual field measurement. 
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Finally, wave data are read into the program and the simulation begins. 
(For information regarding "Read Formats" for the various constants and 
variables, see Appendix E). Wave data required are wave height, wave period, 
wave angle relative to the x-axis of the model at a depth, WDEPTH and the 
duration of the wave climate (HS, T, ALPWIS, and a combination of NTIMES x 
DELT, respectively, in the model). As is always the case with numerical 
models, the time step and space steps are very important to both stability 
and accuracy. Time steps on the order of 3 to 6 hours (10,800 to 21,690 
seconds) or less are recommended. However, the complexity of the bathymetry, 
variation and time series of the wave data, constants used (especially 
Corr) along with several other factors, greatly influences the stability 
and accuracy of the results. 


Table A-1 lists several of the important variables in the computer 
program. 


Table A-1. List of important variables in the program 


ABAND The input banded matrix which stores the values from equation (37) 
ADEAN The value of the scaling parameter in the equilibrium beach profile 
ALPHAS The angle a contour makes with the x-direction base line 


(counter-clockwise is positive) 


ALPWIS The angle (-90° to +90°) the wave crest makes with the 
x-direction (counter-clockwise is positive) 


AMP The amplitude of the diffracted wave in the shadow zone 


ANGGEN The wave angle at a depth, WDEPTH 


ANGLOC The local contour orientation angle 

AWARE See equations (36) and (37) 

BERM The height of the berm above water level 

BMATRX The matrix which, upon solution of the banded matrix problem 
yields the new y values 

C The wave celerity 

CAPPA The breaking wave index 

cc Constant which establishes the width of the distribution of 
sediment transport across the surf zone 

CG The group velocity throughout the wave field 

CGEN The linear wave theory celerity at a depth, WDEPTH 
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CGGEN 
co 
COFF 


CONST 


CONST6 
DEEP 
DEEPB 


DEEPBI 


DELT 
DIAM 
DISTR 
DX 


DY 


HBQ 


The linear wave theory group velocity at a depth, WDEPTH 
The deepwater, linear wave theory wave celerity 


The onshore-offshore transport rate coefficient within the surf 
zone 


The constant in the longshore sediment transport relationship 
(Osu) 


The space step, DX, multiplied by the activity coefficient 
The water depth at any grid location 


The initial breaking depth along each profile (between adjacent 
profiles) 


The initial breaking depth along each profile (at each profile, 
rather than between them) 


The time-step in seconds (DELT x NTIMES = wave condition duration) 
The mean diameter of the sediment particles 
See equations (36) and (37) 


The alongshore space-step in the x-direction (distance between I 
values) 


The onshore-offshore space-step in the y-direction as defined by 
the stepped profile 


The deepwater, linear wave theory wavelength 

The wavelength at the tip of the structure 

The change in the wave number which is acceptably small 
The acceleration of gravity (32.17 feet/second2) 

The specific weight of seawater 

The wave height throughout the wave field 


The maximum wave height which could exist throughout the wave 
field (where H = 0.78 * h) 


The initial breaking wave height along any profile at the y values 
rather than between them 


The initial breaking wave height along any profile, between 
adjacent profiles 
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HGEN Average wave height at a depth, WDEPTH 


HS The significant wave height input 

I The longshore grid location 

IBREAK The leeward side of the initial breaker location J value 

IJET Represents the lesser I value adjacent to the structure (these 
must be evenly spaced alongshore) 

IMAX The total number of grid points in the x-direction (alongshore) 

J The offshore contour location 

JMAX The value of the seawardmost contour simulated 

JUSE (JMAX + 2) the seawardmost contour at which the wave field is 
calculated 

Jl Landward contour of refraction zone 

J2 Seaward contour of refraction zone 

J1REF Landward J values of boundary of refraction zone 

J2REF Seaward J values of boundary of refraction zone 

MMAX The number of shore-perpendicular structures to be simulated 
(present maximum of 16) 

NITER The counterindex in the refraction routine 

NTIME The counterindex in the time simulation "DO" loop 


NTIMES The number of iterations of time-step, DELT, for which a 
particular wave is simulated 


NUNIV The total number of time-steps simulated at any time 

PI The value of 7 = 3.141592654 

POROS The porosity of the sediment 

QX The longshore sediment transport rate at a specific location 

QXTOT The total alongshore sediment transport rate due to the height and 
angle of the initial breaking wave 

QY The onshore-offshore sediment transport rate at a specific location 

R See equations (36) and (37) 
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REPOSE The angle of repose of the sediment 


RHO The mass density of seawater 

RHOND The dimensionless distance from the tip of structure where 
diffraction is initiated 

RHOS The mass density of sediment 

RK The wave number 

$3 See equations (36) and (37) 

SFACE The slope of the shoreface 

SJETTY ie ies of the shore-perpendicular structure (from the base 
ine 

SIGMA The wave radian frequency 

Tf The wave period 

TAU The dissipative interface parameter 

THETA The wave angle throughout the wave field 


THEATO The wave angle at the tip of the structure 


TKSI The longshore sediment transport rate coefficient 

TWOP I Twice the value of 

U See equations (36) and (37) 

UCRIT The critical velocity required to move the sediment according to 
the Sheid's diagram 

V See equations (36) and (37) 

WDEPTH The depth of water in meters to which the input wave conditions 


are to be transformed 


WEQ The equilibrium profile distance between contours as defined by 
the stepped profile 


XCOOR The x-coordinate where the wave field is to be calculated. 
Together with YCOOR, they determine whether the position is within 
or beyond the diffraction shadow zone 

XDISTN The location of the structure along the shoreline in feet 


Y The distance offshore to the contours 
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YCOOR The y-coordinate where the wave field is to be calculated. 
Together with XCOOR, they determine whether the position is within 
or beyond the diffraction shadow zone 


YDISS The value of y after the use of the dissipative interface 
YOLD The previous value of y 

YZERO The berm contour location 

Z1 See equation (37) 

Z2 See equation (37) 
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APPENDIX B 
PROGRAM LISTING 
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100 

200 

300 

400 

500 

600 

700 

800 

900 
1000 
1100 
1200 
1300 
1400 
1500 
1600 
1700 
1800 
1900 
2000 
2100 
2200 
2300 
2400 
2500 
2600 
2700 
2800 
2900 
3000 
3100 
3200 
3300 
3400 
3500 
3600 
3700 
3800 
3900 
4000 
4100 
4200 
4300 
4400 
4500 
4600 
4700 
4800 
4900 
5000 
5100 
5200 
5300 
5400 
5500 
5600 
5700 
5800 
5900 


6000 
6100 
6200 
6300 
6400 
6500 
6600 
6700 
6800 
6900 
7000 
7100 
7200 


C* *** KK KKH PROGRAM IMPLICIT SEDTRAN 

C*THIS PROGRAM IS SET-UP TO HANDLE MULTIPLE GROINS(M<=10). 
COMMON/A/ C(60,20),RK(60,20),Y(60, 20) ,DEEP(60, 20) ,ALPHAS(6O, 20) 
COMMON/AA/YZERO(60) 
COMMON/BB/WEQ(6O0, 20) 
COMMON/B/ THETA(6O,20),QXTOT(60), OLDANG(60,20), DY(60,20) 
COMMON/C/ H(60,20),CG(60, 20) ,HOLD(60, 20) ,HB(60, 20) , YB(60) 
COMMON/N USED/JUSE,T,CO,CGEN, CGGEN, ANGGEN, DX ,BERM, THETAO( 10) , MMAX 
COMMON/D/SIGMA,G,ELO, UMAX, IMAX,PI, TWOPI ,PIO2,HGEN, IJET(10),SUJETTY 
COMMON/F/ADEAN, REPOSE,DIAM 
COMMON/AAA/DELT ,NTIMES 
COMMON/COUNT /NUNIV 
COMMON/EXPL/QYEXP(60, 20), YIMP(60, 20) 
DIMENSION CHANGE(20),HC(10),TC(10) 
DIMENSION YORIG(60, 20), YZEROO(60) ,SANGLE(20) 
NUNIV=O 
JMAX=8 
JUSE=JMAX+2 
IMAX=50 
PI=3.141592654 
TWOPI=PI1*2. 
PIO2=PI/2.0 
REPOSE=32.*TWOPI/360. 
WRITE(6,732) 

732 FORMAT ( TORO OR ORO IO OO OF I OI OI IOI OR a ok rk tk doko 1) 

WRITE(6,733) 

733 FORMAT(2X,’TO WHAT DEPTH ARE THE WAVES TO BE TRANSFORMED’ ) 
C*WDEPTH MUST BE A DEPTH BEYOND THE END OF THE STRUCT, PREFERABLY AT 
C**DEEP(JUMAX) OR GREATER(OR ELSE SNELL’S LAW OR SHOAL COULD BLOWUP IN 
C***DEEPER WATER. IT’S IN METERS HERE! 

READ(5,770) WDEPTH 
770 FORMAT(10X,F 10.3) 
WDEPTH=WDEPTH*3. 28084 
WRITE(6,762) WDEPTH 
762 FORMAT(2X,"THE DEPTH (IN FT) WAVES TRANSFORMED TO, WDEPTH= ". 
* F10.3) 
WRITE(6,732) 
WRITE(6,777) 
777 FORMAT(2X,"ITS TIME FOR SJETTY, BERM, SFACE, AND DIAM",/) 
C*SJETTY MUST BE MUCH LESS THAN Y(I,JUMAX). 
READ(5,776) SUETTY,BERM, SFACE,DIAM 
776 FORMAT(2F10.3,F10.4,F10.3) 
WRITE(6,761) SUETTY 


761 FORMAT(2X,’THE LENGTH OF THE STRUCTURE, SJETTY= ’,F10.3) 
WRITE(6,740) BERM 

740 FORMAT(2X,’THE HEIGHT OF THE BERM, BERM= ‘,F10.3) 
WRITE(6,739)SFACE 

739 FORMAT(2X,’THE SLOPE OF THE BEACH FACE, SFACE= ’,F10.4) 
WRITE(6,738) DIAM 

738 FORMAT(2X,’THE SEDIMENT DIAMETER, DIAM= ’,F10.3) 


WRITE(6,732) 
780 FORMAT(2X,’SUPPLY MMAX( THE NO. OF GROINS) AND THEIR I-LOC’,/) 
UCRIT=16.3*SQRT (DIAM*O.00328) 
C*THE NO. OF MULTIPLE GROINS,MMAX MUST BE GIVEN THEIR X LOCATIONS. 
READ(5,779) MMA X 
779 FORMAT(I3) 
DO 760 M=1,MMAX 
C*IJET REPS LESSER I-VALUE ADJACENT TO STRUCTURE. 
760 READ(5,779) IJET(M) 
WRITE(6,759) (M,IJET(M),M=1,MMAX) 
759 FORMAT(2X,’THE NUMBER’,I5,’ GROIN IS LOCATED AT GRID’,I5) 
WRITE(6,732) 
C*CONVERT TO RADIANS. 
C*FIRST MUST GIVE Y COORS POSITIONS AND DEPTHS. 
C*FIRST, MUST SET UP ALL OF THE DEEP-VALUES. 
WRITE(6,773) 
773 FORMAT(2X,"NOW ENTER THE VALUE OF ADEAN") 
READ(5,774)ADEAN 
774 FORMAT(F10.4) 
WRITE(6,749) ADEAN 
749 FORMAT(2X,’THE VALUE OF ADEAN= ’,F10.4,’ IN THE EQ. H=AY**2/3’) 
WRITE(6,732) 
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7300 WRITE(6,772) 


7400 772 FORMAT(2X,"READ IN THE SPACE STEP,TIMESTEP",/) 

7500 READ(5,775) DX ,DELT 

7600 775 FORMAT(2(F10.3)) 

7700 WRITE(6,737) DX 

7800 737 FORMAT(2X,’THE VALUE OF THE LONGSHORE SPACE-STEP, DX= ’,F10.3) 
7900 WRITE(6,736) DELT 

8000 736 FORMAT(2X,’THE TIME-STEP IN SECONDS, DELT= ’,F10.3) 

8100 DATWARCHANGE/Hm 2 Se Seite dis 6 1405 4725) 39" BOs) 1O*xO.0/, 
8200 DO 220 J=1,UJUMAX+2 

8300 DO 220 I=1,IMAX 

8400 220 DEEP(I,JU)=CHANGE (J) 

8500 DATA(HC(1),1=1,8)/1.87,0.5,0.35, .25, .21,.20,.19,.19/ 

8600 DATAGICCID IM 8))),2os 4p Ga Olin ioeued|4ee/) 

8700 DO 200 J=1,JUMAX+2 

8800 DO 200 I=1, IMAX 

8900 200 Y(1I,U+1)=(0.5*(DEEP(I,JU+1)+DEEP(I,U))/ADEAN) **1.5+Y(I,1) 
9000 WRITE(6,732) 

9100 Co Fe tee OK KOK 

9200 C*WE WILL ALWAYS REQUIRE Y(I,JMAX+2) TO COMPUTE DY AND YBAR. 

9300 C*WE WILL ALWAYS REQUIRE DEEP(I,JMAX+2) TO COMP SEDIMENT TRANSPORT. 
9400 C % Ke ee eK ee eK OK 

9500 WRITE(6,734) 

9600 734 FORMAT(2X,’THE BOUNDARY Y-VALUES, I=1,IMAX ARE AS FOLLOWS’ ,/) 
9700 WRITE(6,801) (Y¥(1,JU),J=1, UMAX+2) 

9800 WRITE(6,801) (Y CIMAX,JU),J=1, UMAX+2) 

9900 WRITE(6,732) 

10000 WRITE(6,735) 

10100 735 FORMAT(/,2X,’THE DEPTHS BETWEEN CONTOURS ARE AS FOLLOWS’ ,/) 
10200 WRITE(6,801) (DEEP(1,JU),J=1, UMAX+2) 

10300 WRITE(6,732) 

10400 801 FORMAT(2X,10(F8.2)) 

10500 DO 2 I=1,IMAX 

10600 2 YZERO(1I)=Y(I,1)-(BERM/SFACE) 


10700 C*WILL COMPUTE THE EQUIL WIDTH BETWEEN CONTOURS, HERE. 
10800 DO 1 I=1,IMAX 


10900 WEQ(I,1)=Y(1,1)-YZERO(I) 

11000 DO 1 J=1,JMAX 

11100 IF(JU.NE.1) GO TO 32 

11200 YTEMP1=0.0 

11300 GO TO 33 

11400 32 YTEMP1=((0.5*(DEEP(I,JU-1)+DEEP(I,U)))/ADEAN)**1.5 
11500 33 YTEMP2=((0.5*(DEEP(1I,J)+DEEP(1I,U+1)))/ADEAN) **1.5 
11600 WEQ(1I,JU+1)=YTEMP2-YTEMP 1 

11700 4 CONTINUE 


11800 C*LET’S STORE THE ORIG VALUES TO COMPUTE VOL CHANGES OVER CONTOURS,LATER 
11900 DO 796 I=1,IMAX+1 


12000 YZEROO(I)=YZERO(I) 

12100 DO 796 J=1,JMAX+2 

12200 796 YORIG(I,J)=Y(I,u) 

12300 CEO OR III I I IOI II III ICICI I aC ICR kok eke kee ae ak a 
12400 C*READ THE DISK FILE AND TRANSFORM PARAMETERS INTO MY UNITS. 

12500 Co ¥ RO RK KR KK | ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! 1 ! ! ! ! ! ! ! ! | eo KK EK KK KK 
12600 C*ALL ADJUSTMENTS TO WAVE ANGLE,HEIGHT,CELERITY,GROUP VEL, WILL BE MADE 
12700 C**HERE, AND THRUOUT THE REST OF THE PROG, THEY WILL BE AS IF OCCURRED 
12800 C***AT WDEPTH! 

12900 798 READ(5,799, END=1000) HS,T,ALPWIS 

13000 799 FORMAT(10X,3F6.1) 

13100 NTIMES=1 

13200 NCHECK=NUNIV+NTIMES 

13300 HGEN=O. 707107 *HS 

13400 SIGMA=TWOPI/T 

13500 G=32.17 

13600 CO=G*T/TWOPI 

13700 ELO=CO*T 

13800 IF(T.LE.2.0) GO TO 797 

43900 HCC=0.23 

14000 DO 444 I=2,7 

14100 T2=TC(I) 

14200 IF(T.GT.T2) GO TO 444 

14300 T1=TC(1-1) 

14400 DELTAT=T2-T1 
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14500 
14600 
14700 
14800 
14900 
15000 
15100 
15200 
15300 
15400 
15500 
15600 
15700 
15800 
15900 
16000 
16100 
16200 
16300 
16400 
16500 
16600 
16700 
16800 
16900 
17000 
17100 
17200 
17300 
17400 
17500 
17600 
17700 
17800 
17900 
18000 
18100 
18200 
18300 
18400 
18500 
18600 
18700 
18800 
18900 
19000 
19100 
19200 
19300 
19400 
19500 
19600 
19700 
19800 
19900 
20000 
20100 
20200 
20300 
20400 
20500 
20600 
20700 
20800 
20900 
21000 
21100 
21200 
21300 
21400 
21500 
21600 


DT=(T-T1)/DELTAT 
DTT=(T2-T)/DELT 
HCC=HC(1)*DT+HC( 1-1) *DTT 
GO TO 446 
444 CONTINUE 
446 CONTINUE 
IF(HGEN.LT.HCC) GO TO 797 
ANGGEN=ALPWIS*TWOPI/360. 
CC FO I I II I OI I OR 2 IK aK oe ok 
CALL WVNUM(WDEPTH, T, DUMKK) 
C*ANGGEN,HGEN,CGEN,CGGEN REPRESENT THE WAVE ANGLE,HEIGHT,CELERITY AND 
C**GROUP VEL(RESPECT.) OF THE SPECIFIED WAVE INPUT AT A DEPTH, WDEPTH 
CALL WVNUM(11.0,T,DUMKKK) 
C11=TWOPI/(T*DUMKKK ) 
CG11=0.5*C11*(1.+(2. *DUMKKK*11.0/SINH(2.*DUMKKK*11.0))) 
CGEN=TWOPI/(T*DUMKK) 
CGGEN=0.5*CGEN*(1.+(2.*DUMKK*WDEPTH/SINH(2.*DUMKK*WDEPTH) ) ) 
CALL TRANS 
797 IF(NCHECK.NE.NUNIV) NUNIV=NCHECK 
709 GO TO 798 
1000 CONTINUE 
STOP 
END 
Co KR ee Eo oo oo a oa oo eo oo I oo KK 2 KK oR RK eK ok Kk OK 
SUBROUTINE TRANS 
C*THIS SUBROUTINE WILL COMPUTE SEDIMENT TRANSPORT 
COMMON/A/ C(60,2C),RK(60,20),Y(60, 20) ,DEEP(60, 20) , ALPHAS(60, 20) 
COMMON/AA/YZERO(60) 
COMMON/BB/WEQ(60, 20) 
COMMON/B/ THETA(60,20),QXTOT(60), OLDANG(60,20), DY(60,20) 
COMMON/C/ H(60,20),CG(60, 20) ,HOLD(60, 20) ,HB(60, 20), YB(60) 
COMMON/N USED/JUSE,T,CO,CGEN, CGGEN, ANGGEN, DX ,BERM, THETAO( 10) , MMAX 
COMMON/D/SIGMA,G,ELO, JMAX, IMAX,PI, TWOPI ,PIO2,HGEN, IJET( 10), SUETTY 
COMMON/E/RHO,RHOS,POROS,CONST,TKSI 
COMMON/F/ADEAN, REPOSE,DIAM 
COMMON/G/IBREAK(60) , HNONBR( 20) 
COMMON/P/HBQ(60) , DEEPB(60) 
COMMON/ZZZ/NTIME 
COMMON/AAA/DELT,NTIMES 
COMMON/COUNT/NUNIV 
DIMENSION YOLD(60, 20),R(60,20),S(60, 20) ,HC(60, 20) ,QY(60, 20), YDISS( 
* 60,20) 
DIMENSION RHS1(60, 20) ,S3(60, 20), THETAB(60, 20) , ANGLOC(60, 20) 
DIMENSION DISTR(60, 20) ,AWARE(60, 20 


CE OK eo ko I a I I Ra IC ko ke ake ke ok ok FIC ok ok oo oko akc ok ke ake ok ok ke ke ok ok ke ke a ke ake ake ok ake ake ke ok ake 
Co ee eo ke ko i ok ke oko ok 2 oo ie ke ke kk kok oko EE oe ke keke ke ake ke ke ok ke ok oe ok oe oe ok oe ok ok ote ote ote ote oe oe ake ake ake ake ke ke ke 


CR eee ee ee ee NO ITE: : SIZE OF ABAND AND XL HAVE TO BE CHANGED 
CRE OR ERE eee Le ACCORDING TO JMAX+1+JMAX AND JMAX+1,RESPECT. 
CORSE EAE AE RE SEICOE AE REAR BERS CHANGE REQ’D AT 7040 AND 18650 


CK I ok OK I 2K a aK ook koko ook ok ok ok oka oko kok oe oe ok ok ok oe ok ok ok ok ke kok ke ke oe ke 


* ), BMATRX (432), ABAND(432, 19) ,QX(60, 20) ,XL(432, 10) , CONST6(60, 20) 
COMMON/MP/ RKB(60),HBI(60) ,DEEPBI(60) 
COMMON/EXPL/QYEXP(60, 20), YIMP(60, 20) 

DIMENSION SANGLE(20) 
C*LET’S ZERO-OUT ALL OF THE DIMENSIONED MATRICES. 

DO 1000 J=1,UMAX+2 

SANGLE(JU)=0.0 

DO 1000 I=1,IMAX+2 

YOLD(I,JU)=0.0 

R(I,JU)=0.0 

S(I,JU)=0.0 

HC(1I,JU)=0.0 

QY(I,JU)=0.0 

YDISS(I,U)=0 

RHS1(1I,U)=0. 

S$3(1,JU)=0.0 
THETAB(I,U) 
ANGLOC(I,Ju) 
DISTR(I, u)= 
AWARE(I,J)= 
QX(I,JU)=0.0 
CONST6(I,JU)=0.0 


O. 
oO. 
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21700 
21800 
21900 
22000 
22100 
22200 
22300 
22400 
22500 
22600 
22700 
22800 
22900 
23000 
23100 
23200 
23300 
23400 
23500 
23600 
23700 
23800 
23900 
24000 
24100 
24200 
24300 
24400 
24500 
24600 
24700 
24800 
24900 
25000 
25100 
25200 
25300 
25400 
25500 
25600 
25700 
25800 
25900 
26000 
26100 
26200 
26300 
26400 
26500 
26600 
26700 
26800 
26900 
27000 
27100 
27200 
27300 
27400 
27500 
27600 
27700 
27800 
27900 
28000 
28100 
28200 
28300 
28400 
28500 
28600 
28700 
28800 
28900 


4000 CONTINUE 
RHO=1.99 
RHOS=5. 14 
POROS=0.40 
CONST=0.77 
CAPPA=0.78 
TAU=0.25 
TKSI=(CONST*RHO*SQRT(G) )/((RHOS-RHO)*(1.O-POROS)*16.0*SQRT(CAPPA) 
C* QX(I,JU) IS THE TRANSPORT BETWEEN THE (1I,JU+1) AND (I,J) CONTOURS. 
C*THE ‘DO 1 LOOP’ SIMULATES TIME---TIME=DELT*NTIMES. 
COFF=0.00001 
GAMMA=RHO*G 
DO 1 NTIME=1,NTIMES 
NUNIV=NUNIV+1 
C*THE MATRICES ABAND AND BMATRX MUST BE "ZEROED OUT" 
K=O 
DO 26 I=2,IMAX-1 
DO 26 J=i,JUMAX 
K=K+1 
BMATRX(K)=0.0 
DO 26 L=1,JUMAX+1+JMAX 
26 ABAND(K,L)=0.0 
XNT IME=1.0*(NTIME) 
CALL PREDIF 
C*SMOOTHING OF THE WAVE ANGLE,THETA, IS RE’D TO ACCT FOR DIFF EFFECTS. 
CALL SMOOTH( THETA, IMAX, UMAX, IJET, SUETTY ,MMAX,Y) 
CALL QTRAN 
C*FIRST THE LONGSHORE SEDIMENT TRANSPORT WILL BE DISTRIBUTED 
C****ACROSS THE SURF ZONE.... 
CC=Hh25 
C***QX(I,U) WILL BE DETERMINED BY SUBTRACTING FROM THE INTEGRAL 
C**OF QX FROM DEEP(I,J-1) TO INFINITY, THE INTEGRAL OF QX FROM DEEP(I,u) 
C***TO INFINITY. IN THIS WAY THE SEDIMENT TRANS FROM JMAX OUT GETS 
C***INCLUDED IN QX(I,JMAX). TO INCLUDE THE SWASH TRANS, WHEN J=1 
C*WE WILL SUBTRACT FROM 2 TO INFINITY FROM 1.0 
C*LOOP FOR VALUES WHICH ARE HELD CONST AND STORED. 
THETAB(1,1)=0.5*(THETA(1,1)+0.0) 
R(1,1)=0.5/(DX*(DEEP(1,1)+BERM/2.)) 
DO 290 I=2,IMAX 
R(1,1)=0.5/(DX*(DEEP(1,1)+BERM/2.)) 
c* THETAB(I,1)=0.25*(THETA(I,1)+THETA(I-1,1)+0.+0. ) 
THETAB(1I,1)=0.5*(THETA(I,1)+THETA(I-1,1)) 
C*NO NEED TO COMPUTE PROP ANGLE AT STRUCTS BECAUSE QX =0.0O AT IJET(M)+1 
ANGLOC(I,1)=ATAN((Y(1,1)-Y(1I-1,1))/DX) 
C*HBQ(IJET(M)+1) IS PROPERLY SET IN THE SUBROUTINE QTRAN. 
DISTR(I,1)=1.0-EXP(-((DEEP(I, 1)**1.5+HBQ(1I)*ADEAN**1.5)/ 
* (CC*DEEPB(1)**1.5))**3) 
DISTR(I,1)=DISTR(I,1)*TKSI*HBOQ(I)**2.5 
DO 280 J=2,JUMAX 
R(1,JU)=0.5/(DX*(DEEP(1I,U)-DEEP(I,U-1))) 
THETAB(I,JU)=O.5*(THETA(I,JU)+THETA(I-1,uU)) 
ANGLOC(I,J)=ATAN((Y(1I,JU)-Y(1I-1,U))/DX) 
DISTR(I,JU)=EXP(-((DEEP(I,JU-1)**1.5+HBQ(1I) *ADEAN**1.5)/(CC*DEEPB(1) 
* **4.5))**3)-EXP(-((DEEP(1I,JU)**1.5+HBQ(I) *ADEAN**1.5)/(CC* 
* DEEPB(I)**1.5))**3) 
DISTR(I,JU)=DISTR(I,JU)*TKSI*HBQ(1)**2.5 
290 CONTINUE 
DO 301 J=1,JMAX 
DO 301 I=2,IMAX 
AWARE(I,JU)=DELT*R(I,JU)*(QX(I,JU)-QX(I+1,JU)+QY(1I,JU)-QY(I,U+1))+Y(I,U 
* 
) 
S1=2.*SIN(THETAB(I,JU))*COS(THETAB(I,JU))*(-1.+2.*(COS( 
* ANGLOC(I,uJ)))**2) 
S2=COS(2.*THETAB(I,U))*COS(ANGLOC(1I,U))/(SQRT(DX**2+ 
* (Y(I,U)-Y(1-1,J))**2)) 
$3(1,JU)=S2*DISTR(I,u) 
IF(SUJETTY.EQ.0.0) GO TO 302 
DO 325 M=1,MMAX 
IF(I.NE.IJET(M)+1) GO TO 325 
IF(THETAO(M).GE.O.0) ISIDE=IJET(M) 
IF(THETAO(M).LT.O.0) ISIDE=IJET(M)+1 
YSEA=0.5*(Y(ISIDE,J)+Y(ISIDE,uU+1)) 
YSHORE=0.5*(Y(ISIDE,JU)+Y(ISIDE,uU-1)) 
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29000 
29100 
29200 
29300 
29400 
29500 
29600 
29700 
29800 
29900 
30000 
30100 
30200 
30300 
30400 
30500 
30600 
30700 
30800 
30900 
31000 
31100 
31200 
31300 
31400 
31500 
31600 
31700 
31800 
31900 
32000 
32100 
32200 
32300 
32400 
32500 
32600 
32700 
32800 
32900 
33000 
33100 
33200 
33300 
33400 
33500 
33600 
33700 
33800 
33900 
34000 
34100 
34200 
34300 
34400 
34500 
34600 
34700 
34800 
34300 
35000 
35100 
35200 
35300 
35400 
35500 
35600 
35700 
35800 
359300 
36000 
36100 
36200 


IF(YSEA.GT.SJETTY.AND.YSHORE.GT.SJETTY) GO TO 302 
IF(YSEA.GT.SJETTY.AND.YSHORE.LE.SJETTY) GO TO 298 
C*BECAUSE A NO FLOW B.C. IS USED ALONG THE STRUCT, NO ATTN WAS PAID 
C**TO GETTING PROPER VALUES OF ANGLOC, THETAB,DISTR,ETC. 
S$3(1,JU)=0.0 
DISTR(I,JU)=0.0 
GO TO 302 
325 CONTINUE 
GO TO 302 
C***ABOVE, ALL PARAMETERS(I.E.,S1,S2,S3, THETAB,DISTR,ANGLOC) 
C***ARE COMPUTED AS IF THE STRUCT IS NOT THERE. THE B.C. AT THE 
C***STRUCT TIP ASSUMES QX COMPUTED AS IF NO STRUCT PRESENT AND THEN 
C***BYPASSES ACCORDING TO "RATIO". 
298 RATIO=(YSEA-SJETTY)/(YSEA-YSHORE ) 
S$3(1,JU)=S3(1I,JU)*RATIO 
DISTR(1I,JU)=DISTR(I,U)*RATIO 
302 RHS1(I,U)=DISTR(I,JU)*S1-S3(1,U)*(Y(1,U)-Y(1-1,u)) 
301 CONTINUE 
CALL BREAK(IMAX, JUMAX) 
C*TO DETERMINE DECAY OF CONST6(I,JU),NEED WAVE NO. AT BREAKING. 
DO 754 I=1,IMAX+1 
754 CALL WVNUM(DEEPBI(I),T,RKB(I)) 
C*USING SHIELD’S DIAG,Y AXIS=0.05 & (TAUO=RHO*C*U**2),GET UCRIT(FT/SEC) 
UCRIT=16.3*SQRT(DIAM* .00328 ) 
DO 750 I=1, IMAX+14 
CONST6(I, 1)=COFF*DX 
DO 750 J=2,UMAX+2 
C*CONST6(1I,U) GOES W/ QY(I,JU) WHICH IS ASSOC W/ DEEP(I,uU-1) 
IF(DEEP(I,J-1).LE.DEEPBI(I)) GO TO 751 
C*HERE, MUST CAUSE COFF TO DECAY (WE’RE BEYOND SURF ZONE) 
UMAXB=HBI (1 )*G*T*RKB(1I)/(2.*TWOPI*COSH(RKB(1I)*DEEPBI(I))) 
UMAX=H( I, JU-1)*G*T*RK(I,J-1)/(2. *TWOPI*COSH(RK(I,U-1)*DEEP(I,U-1))) 
IF(UCRIT.LT.UMAX.AND.UCRIT.LT.UMAXB) GO TO 749 
CONST6(I,U)=0.0 
GO TO 750 
749 TOP=0.01*H(1I,JU-1)**3*SIGMA**3/(SINH(RK(I,J-1)*DEEP(I,JU-1))**3) 
BOT=DEEP(I,U-1)*(0.625*TWOPI*G**1.5*O.78**2*ADEAN**1.5+ 
*(0.01*HBI(1)**3*SIGMA**3/(DEEPBI(1I)*(SINH(RKB(I)*DEEPBI(1)))**3))) 
CONST6(1,JU)=DX*COFF*TOP/BOT 
GO TO 750 
751 CONST6(1I,U)=COFF*DX 
750 CONTINUE 
K=O 
C**PUT INTO BANDED FORM USING THE ALGORITHM A(M,N)->B(M,NN) WHERE 
C***NN=KB+1-M+N(KB IS THE NUMBER OF LOWER CODIAGONALS(=UMAX,HERE)). 
DO 304 I=2,IMAX-1 
DO 304 J=1,UMAX 
K=K+1 
AWARE(1,U)=AWARE(I,JU)+DELT*RHS1(1,JU)*R(I,U)-DELT*R(I,JU)*RHS1(I+1,J 
as )+DELT*R(I,U)*CONST6(1I,JU)*WEQ(1,JU)-DELT*R(I,JU)*CONST6(I,JU+1)* 
* WEQ(I,U+1) 
YDUM=YZERO(1) 
IF(JU.NE. 1) YDUM=Y(I,J-1) 
AWARE(I,U)=AWARE(1I,JU)+DELT*R(I,JU)*CONST6(I,JU)*O0.5*(YDUM-Y(I,U)) 
* -DELT*R(I,U) *CONST6(I,JU+1)*O0.5*(Y(1I,JU)-Y(1I,U+1)) 
U=DELT*R(I,JU)*S3(I,u) 
V=DELT*R(I,JU)*S3(I+1,U) 
Z1=DELT*R(I,JU)*CONST6(I,JU)*0.5 
Z2=DELT*R(I,JU)*CONST6(I,U+1)*O.5 
C*NOW WILL SET UP THE MATRICES ABAND AND BMATRX. 
ABAND(K, UMAX+1)=1.O0+U+V+Z1+Z2 
IF(I.NE.2) GO TO 305 
AWARE(1I,U)=AWARE(I,JU)+U*Y(I-1,U) 
GO TO 310 
305 ABAND(K,1)=-U 
310 IF(I.NE.IMAX-1) GO TO 306 
AWARE(1,JU)=AWARE(I,JU)+V*Y(IMAX,J) 
GO TO 311 
306 ABAND(K, JMAX+1+JMAX )=-V 
311 IF(JU.NE.1) GO TO 307 
ABAND(K, JMAX+1)=ABAND(K, UMAX+1)-Z1 
AWARE(1,1)=AWARE(I,1)+Z1*(YZERO(I)-Y(I,1)) 
GO TO 312 
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36300 307 ABAND(K, JUMAX)=-Z1 


36400 312 IF(U.NE.UMAX) GO TO 308 

36500 AWARE(1,JU)=AWARE(1I,JU)+Z2*Y(I, UMAX+1) 

36600 GO TO 309 

36700 308 ABAND(K,JMAX+2)=-Z2 

36800 309 BMATRX(K)=AWARE(I,J) 

36900 304 CONTINUE 

37000 KMAX=K 

37100 C**CALL IMSL ROUTINE LEQT1B TO SOLVE THE BANDED MATRIX. 

37200 CALL LEQT1B(ABAND,KMAX, JUMAX , UMAX,432,BMATRX,1,432,0,XL, IER) 
37300 C*NOW, GIVE Y’S THEIR NEW VALUES STORING OLD VALUES IN YOLD. 
37400 K=O 

37500 DO 315 I=2,IMAX-1 

37600 YOLD(1I, UMAX+1)=Y(1I, UMAX+1) 

37700 DO 315 J= 1,JMAX 

37800 K=K+1 

37300 YOLD(I,J)=Y(I,U) 

38000 Y(1I,JU)=BMATRX(K) 

38100 315 CONTINUE 

38200 DO 320 J=1,UMAX+3 

38300 YOLD(1,U)=Y(1,u) 

38400 320 YOLD( IMAX,JU)=Y(IMAX,J) 

38500 C*WILL USE ABBOTT’S DISSIPATIVE INTERFACE TO RID HIGH FREQ OSCILLATIONS 
38600 DO 650 J=1,JUMAX 

38700 DO 650 I=2,IMAX-1 

38800 YDISS(1I,JU)=TAU*Y(I-1,JU)+(1.-2.*TAU)*Y(I,JU)+TAU*Y(I+1,U) 
38900 IF(SJETTY.EQ.0.0) GO TO 650 

39000 DO 649 M=1,MMAX 

39100 IF(I.NE.IJET(M).AND.I.NE.IJET(M)+1) GO TO 649 

39200 IF(Y(IJVET(M),JU).GT.SJUETTY.OR.Y(IJUET(M)+1,U).GT.SJETTY)GO TO 649 
39300 IF(I.EQ.IJET(M) )YDISS(I,JU)=TAU*Y(I-1,U)+(1.-TAU)*Y(I,J) 
39400 IF(1I.EQ. (IUET(M)+1))YDISS(1I,U)=TAU*Y(I+1,JU)+(1.-TAU)*Y(I,U) 
39500 649 CONTINUE 

39600 650 CONTINUE 

39700 DO 651 J=1,JUMAX 

39800 DO 651 I=2,IMAX-14 

39900 651 Y(I,JU)=YDISS(I,u) 


40000 C*THIS LOOP WILL STORE THE IMPLICIT Y VALUES REQ’D TO COMP QY&QX 
40100 DO 40 I=1,1IMAX+1 


40200 DO 40 J=1,JUMAX+3 

40300 40 YIMP(I,J)=Y(I,U) 

40400 C*THIS LOOP WILL EXPLICITLY MOVE CONTOURS SEAWARD IF REPOSE EXCEEDED. 
40500 KOUNT =O 

40600 SLOPEM=TAN(O.9*REPOSE ) 

40700 DO 48 I=1,IMAX 

40800 43 KOUNT=KOUNT+1 

40300 IF (KOUNT . GT .50000) GO TO 41 


41000 C*LET US COMPUTE ALL THE SLOPES(PSLOP) FOR EACH CHANGE IN DEPTH. 
41100 DO 47 J=1,JMAX+1 


41200 DUM=-BERM/2.0 

41300 IF(JU.NE.1) DUM=DEEP(I,U-1) 

41400 DELH=0.5*(DEEP(I,JU+1)+DEEP(I,JU))-0.5*(DEEP(I,U)+DUM) 

41500 PSLOP=DELH/(Y(I,JU+1)-Y(1I,U)) 

41600 47 SANGLE(J)=ATAN(PSLOP ) 

41700 C*FIND THE MIN NEG SLOPE ANGLE OR THEN THE POS SLOPE>REPOSE OR FORGET IT 
41800 ASLOPM=-1.0E50 

413900 ASLOPP=0.0 

42000 DO 46 J=1,JMAX+1 

42100 IF(SANGLE(JU).GT.0.0) GO TO 45 

42200 IF(SANGLE(JU).GT.ASLOPM)ASLOPM=SANGLE (J) 

42300 TF (ASLOPM.EQ.SANGLE(u)) JM=J 

42400 GO TO 46 

42500 45  IF(SANGLE(J).GT.REPOSE.AND.SANGLE(J).GT.ASLOPP )ASLOPP=SANGLE (J) 
42600 IF(ASLOPP.EQ.SANGLE(Uu)) JP=J 

42700 46 CONTINUE 

42800 IF(ASLOPM.EQ.-1.0£50.AND.ASLOPP.EQ.0.0) GO TO 42 

42900 IF(ASLOPM.EQ.-1.0E50) GO TO 44 

43000 DUM=-BERM/2. 

43100 IF(JUM.NE. 1) DUM=DEEP(I,UM-1) 

43200 ALTER=((0.5/SLOPEM*(DEEP(1I,JM+1)-DUM) )-(Y(I,JUM+1)-Y(I,UM)))/ 
43300 * (1.0+((DEEP(I,JM+1)-DEEP(I,UM))/(DEEP(I,JUM)-DUM) ) ) 

43400 Y(1I,JUM+1)=Y(1,UM+1)+ALTER 


ll 


43500 
43600 
43700 
43800 
43900 
44000 
44100 
44200 
44300 
44400 
44500 
44600 
44700 
44800 
44900 
45000 
45100 
45200 
45300 
45400 
45500 
45600 
45700 
45800 
45900 
46000 
46100 
46200 
46300 
46400 
46500 
46600 
46700 
46800 
46900 
47000 
47100 
47200 
47300 
47400 
47500 
47600 
47700 
47800 
47900 
48000 
48100 
48200 
48300 
48400 
48500 
48600 
48700 
48800 
48900 
49000 
49100 
49200 
49300 
49400 
49500 
49600 
49700 
49800 
49900 
50000 
50100 
50200 
50300 
50400 
50500 
50600 
50700 


44 


42 
48 


Y(1,JUM)=Y(1,JUM)-(ALTER*(DEEP(I,JUM+1)-DEEP(1I,JUM))/(DEEP(I,UM)-DUM) ) 
QYEXP(1,JM+1)=QYEXP(I,JUM+1)+DX/DELT*ALTER*(DEEP(I , UM+1)-DEEP(I,UM) 
* 
) 

GO TO 43 

CONTINUE 

DUM=-BERM/2. 

IF(JUP.NE.1) DUM=DEEP(I,UP-1) 
ALTER=((0.5/SLOPEM*(DEEP(I,JUP+1)-DUM) )-(Y(I,JUP+1)-Y(I,UP)))/ 

* (1.0+((DEEP(1I,JP+1)-DEEP(1,JUP))/(DEEP(1,UP)-DUM) )) 
Y(I,JUP+1)=Y(1,JUP+1)+ALTER j 
Y(I,UP)=Y(1I,JUP)-(ALTER*(DEEP(1I,JP+1)-DEEP(I,UP))/(DEEP(1I,JUP)-DUM) ) 
QYEXP(1I,UP+1)=QYEXP(1I,JUP+1)+DX/DELT*ALTER*(DEEP(I,UP+1)-DEEP(1,UP) 

* 

) 

GO TO 43 

WEQ(I, UMAX+1)=Y (I, JMAX+1)-Y(1I, UMAX) 

CONTINUE 


C*IF WE GET SENT HERE, LOOP 444 WILL CATCH THE CROSSED CONTOURS. 


319 


318 


323 


CONT INUE 
WE CAN COMPUTE QX’S AND QY’S! 

DO 318 I=2,IMAX 

IMPLIC AND EXPLIC MOVEMENT OF YZERO WILL BE TAKEN CARE OF HERE 
QY(1, 1)=-BERM*DX*(Y(1I,1)-YOLD(I,1))/DELT 
YZERO(I)=YZERO(1I)+(Y(1,1)-YOLD(I,1)) 

DO 318 J=1,JUMAX 
QX(1,U)=RHS1(1,U)-S3(1,JU)*YIMP(I,JU)+S3(1I,U)*YIMP(I-1,u) 
QY(1,U+1)=CONST6(I,JU+1)*(0.5*(YIMP(I,JU)+YOLD(I,JU)-YIMP(I,J+1) 
* -YOLD(I,uU+1))+WEQ(I,uU+1)) 

DO 323 J=1,JUMAX 

QX(1,U)=QX(2,U) 

QX (IMAX+1,JU)=QX(IMAX, J) 


C*TOTAL QYS WILL BE COMP FROM IMPLIC AND EXPLIC VALUES.THEN ZERO QYEXP 


39 


DO 39 I=1,IMAX+1 
DO 39 J=1,UMAX+3 
QY(I,JU)=OY(1I,JU)+QYEXP(I,J) 
QYEXP(I,JU)=0.0 


C*THIS CHECK WILL BOMB THINGS OUT IF CONTOURS HAVE CROSSED. 


DO 444 II=1, IMAX 
DO 444 JJ=1,JUMAX 


C*IF CONTOURS CROSS AT ANY TIME WANT PROGRAM TO STOP! 


150 


151 


152 
103 


19 


444 


C*THE 


C*LET 


926 
800 
Cc* 
c*900 
Cc* 
Cc*903 
Cc* 
C*906 
c* 
Cx755 
801 


107 


15 


IF(Y(II,JUJU).LT.Y(II,JJU+1)) GO TO 444 
WRITE(6, 103) 
WRITE(6, */) NUNIV 
DO 150 J=1,JUMAX 
WRITE(6, 100) (QX(I,JU),1=1, IMAX) 
DO 151 J=1,JUMAX 
WRITE(6, 101) (QY(I,U),1=1, IMAX) 
DO 152 J=1,JMAX 
WRITE(6, 100) (Y(1,JU),1=1, IMAX) 
FORMAT(2X,’THE CONTOURS HAVE CROSSED AND SOMETHING IS WRONG’ ,/) 
DO 19 J=1,UMAX 
WRITE(6,100) (YOLD(I,JU),1=1,IMAX) 
GO TO 445 
CONT INUE 
WRITE(6,*/) NUNIV 
FOLLOWING STATEMENT DETERMINES AT WHAT FREQ EVERYTHING IS WRITTEN! 
IF (MOD(NUNIV, 10) .NE.O) GO TO 1 
’S WRITE ALL OF IT OUT. 
WRITE(6,926) NUNIV 
FORMAT(2X,’THE TOTAL ELAPSED NUMBER OF TIME-STEPS. NUNIV= ’,15,/) 
FORMAT(2X,14(F8.4)) 
DO 900 I=1,IMAX 
WRITE(6,800) (THETA(I,JU),J=1, UMAX) 
DO 903 J=1,UMAX+1 
WRITE(6,801) DEEP(1,uU) 
DO 906 I=1,IMAX 
WRITE(6,800) (H(I,JU),J=1,JUMAX) 
DO 755 J=1,JUMAX 
WRITE(6,800) (CONST6(I,U),1I=1, IMAX) 
FORMAT(2X,14(F8.2)) 
WRITE(6, 107) 
FORMAT(/,2X,’THE LONGSHORE TRANSPORTS,QX, FOLLOW’ ) 
DO 15 J=1,JMAX 
WRITE(6, 100) (QX(1,JU),1=1, IMAX) 


UZ 


50800 
50900 
51000 
51100 
51200 
51300 
51400 
51500 
51600 
51700 
51800 
51900 
52000 
52100 
52200 
52300 
52400 
52500 
52600 
52700 
52800 
52900 
53000 
53100 
53200 
53300 
53400 
53500 
53600 
53700 
53800 
53900 
54000 
54100 
54200 
54300 
54400 
54500 
54600 
54700 
54800 
54900 
55000 
55100 
55200 
55300 
55400 
55500 
55600 
55700 
55800 
55900 
56000 
56100 
56200 
56300 
56400 
56500 
56600 
56700 
56800 
56900 
57000 
57100 
57200 
57300 
57400 
57500 
57600 
57700 
57800 
57900 
58000 


WRITE(6, 108) 
108 FORMAT(/,2X,’THE ON-OFFSHORE TRANSPORTS, QY, FOLLOW’ ) 
DO 17 J=1,JUMAX 
17 WRITE(6, 101) (QY(I,JU),1=1, IMAX) 
WRITE(6, 109) 
109 FORMAT(/,2X,’THE NEW CONTOUR VALUES, Y, FOLLOW’) 
DO 18 J=1,JMAX 
18 WRITE(6, 100) (Y(1I,JU),1=1, IMAX) 
100 FORMAT(2X,13(F9.3)) 
101 FORMAT(2X,13(F9.4)) 
1 CONT INUE 
RETURN 
GO TO 446 
445 STOP 
446 CONTINUE 
END 
CK ee oo eee kk ok OK ka KK KK kK ok ak aK ok kk a ok ie akc ok ok ke ok G 
SUBROUTINE QTRAN 
C*THIS SUBROUTINE CALCS THE BREAKER HEIGHT FOR EACH 
C*OF THE I GRID LINES. METHOD--FINDS Y-LOCATIONS BEFORE AND AFTER 
C*BREAKING HAS OCCURRED BY ‘’REFRAC’, THEN USES SHOALING TO GET THE 
C*HBQ.SNELL’S LAW IS USED FOR REFRACTION OVER THE SHORT DIST TO BREAKING 
C* QX(I,Ju) IS THE TRANS BETWEEN(I-1,JU) AND (I,J) AT THE BLOCKCENT 
COMMON/A/ C(60,20),RK(60,20),Y(60, 20) ,DEEP(60, 20) , ALPHAS(60, 20) 
COMMON/AA/YZERO(60) 
COMMON/B/ THETA(60,20),QXTOT(60), OLDANG(60,20), DY(60, 20) 
COMMON/C/ H(60,20),CG(60, 20) ,HOLD(60, 20) ,HB(60, 20), YB(60) 
COMMON/N USED/JUSE,T,CO,CGEN,CGGEN, ANGGEN,DX,BERM, THETAO( 10) , MMAX 
COMMON/D/SIGMA,G,ELO, JUMAX, IMAX,PI, TWOPI,PIO2,HGEN, IJET( 10), SUETTY 
COMMON/G/ IBREAK(60) , HNONBR( 20) 
COMMON/E/RHO, RHOS,POROS, CONST, TKSI 
COMMON/P/HBQ(60) ,DEEPB(60) 
CAPPA=0.78 
DO 1 I=2,IMAX 
DO 2 JJ=1,JMAX 
JU=JMAX-Ju+1 
HDUM=(H(1I,JU)+H(I-1,JU))*0.5 
HBDUM=(HB(1I,JU)+HB(I-1,JU))*0.5 
C*CAN ONLY USE COND ON ONE SIDE OF STRUCT. CAN’T AVG HERE! 
IF(SJETTY.EQ.0.0) GO TO 4 
DO 4 M=1,MMAX 
IF(I.NE.IVET(M)+1) GO TO 4 
IF(THETAO(M).GE.0.0) ISIDE=IJET(M) 
IF(THETAO(M) .LT.O.0) ISIDE=IJET(M)+1 
C***B.C. AT STRUCT TIP ASSUMES GX COMP AS IF NO STRUCT IS PRESENT. 
YSEA=0.5*(Y(ISIDE,JU)+Y(ISIDE,uJ+1)) 
IF(YSEA.GT.SJETTY) GO TO 3 
HDUM=H(ISIDE,u) 
HBDUM=HB(ISIDE,uU) 
GO TO 3 
4 CONTINUE 
3. IF(HDUM.€T.HBDUM) GO TO 2 
DEEPB(I)=((0.5*(H(I,U+1)+H(I-1,U+1)))*((0O.5*(DEEP(1I,u+1) 
* +DEEP(I-1,JU+1)))**O.25)/CAPPA)**O.8 
HBQ(I)=CAPPA*DEEPB(I) i 
C*HBQ(I) AND DEEPB(I) WILL BE COMPUTED ACCORDING TO THE WAVE DIR. 
C** AT THE STRUCTURE TIP,THETAO. 
IF(SUETTY.EQ.0.0) GO TO 1 
DO 6 M=1,MMAX 
IF(I.NE.IJET(M)+1) GO TO 6 
C**THE TRANSPORTING WAVES WILL BE COMPUTED USING THE WAVE TO PROP SIDE. 
IF(THETAO(M).GE.0O.0) GO TO 11 
DEEPB(1)=(H(IJET(M)+1,uU+1)*DEEP(IJET(M)+1,U+1)**O.25/CAPPA)**0.8 
IBREAK(I)=IBREAK(IJET(M)+1) 
GO TO 12 
41 DEEPB(I)=(H(IVET(M) ,J+1)*DEEP(IJET(M),J+1)**O.25/CAPPA)**0.8 
IBREAK(1I)=IBREAK(IJET(M) ) 
12 HBQ(1)=DEEPB(I)*CAPPA 
GO TO 1 
6 CONTINUE 
GO TO 1 
CONT INUE 
1 CONT INUE 
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58100 
58200 
58300 
58400 
58500 
58600 
58700 
58800 
58900 
59000 
59100 
59200 
59300 
59400 
59500 
59600 
59700 
59800 
59900 
60000 
60100 
60200 
60300 
60400 
60500 
60600 
60700 
60800 
60900 
6 1000 
61100 
61200 
61300 
61400 
61500 
61600 
61700 
61800 
61900 
62000 
62100 
62200 
62300 
62400 
62500 
62600 
62700 
62800 
62900 
63000 
63100 
63200 
63300 
63400 
63500 
63600 
63700 
63800 
63900 
64000 
64100 
64200 
64300 
64400 
64500 
64600 
64700 
64800 
64900 
65000 
65100 
65200 


C*IF THE OFFSHORE WAVE HT IS ZERO, NEVER GET TO HERE. 
C*HOWEVER IF THE H IS SUCH THAT IT WOULD BREAK INSHORE OF Y(I,2) 
C*DEEPB(I) WOULD STILL BE ZERO AND DISTR(I,JU) WOULD BLOW-UP. 
DO 20 I=1,IMAX 
IF (DEEPB(I).GT.0O.0O) GO TO 20 
DEEPB(I)=(H(1I,1)*DEEP(1,1)**0.25/CAPPA)**0.8 
HBQ( I )=CAPPA*DEEPB(1) 
20 CONTINUE 
HBQ(1)=HBQ(2) 
HBQ( IMAX+1)=HBQ( IMAX) 
DEEPB(1)=DEEPB(2) 
DEEPB( IMAX+1)=DEEPB(IMAX) 
RETURN 
END 
CF kt tee ee eK keke te te eke kk a ke ok ke eke kek ke eo oe kk ao ee ek 
SUBROUTINE BREAK(IMAX,UJMAX) 
C*ROUTINE WILL DETERMINE HB AND DEEPB ON THE GRID LINES RATHER 
C* THAN BETWEEN THEM. REQ’D FOR COFF BEYOND SURF ZONE. 
COMMON/A/ C(60,20),RK(60,20),Y(60, 20) ,DEEP(60, 20) , ALPHAS(60, 20) 
COMMON/C/ H(60,20),CG(60, 20) ,HOLD(60, 20) ,HB(60, 20), YB(60) 
COMMON/MP/ RKB(60),HBI(60) ,DEEPBI(60) 
CAPPA=0.78 
DO 1 I=2,IMAX 
DO 2 JJ=1,UJUMAX 
J=JMAX-Ju+1 
IF(H(1I,J).LT.HB(I,U)) GO TO 2 
DEEPBI(1I)=((H(1,U+1)*DEEP(I,JU+1)**O.25)/CAPPA)**0.8 
HBI(I)=CAPPA*DEEPBI(1) 
C***ONCE THE HEIGHT & DEPTH AT BREAKING ARE FOUND, GO TO NEXT GRID-LINE. 
GO TO 1 
2 CONTINUE 
1 CONTINUE 
DO 20 I=1,IMAX 
IF(DEEPBI(I).GT.0.0) GO TO 20 
DEEPBI(1)=(H(1I,1)*DEEP(1,1)**O.25/CAPPA)**0.8 
HBI(1I)=CAPPA*DEEPBI(I) 
20 CONTINUE 
DEEPBI(1)=DEEPBI(2) 
DEEPBI( IMAX+1)=DEEPBI(IMAX) 
HBI(1)=HBI(2) 
HBI ( IMAX+1)=HBI(IMAX) 
RETURN 
END 
CORO OR ORO OR III OR I I IO IOI ICI ICC IIE IIE A IIR Ik a kk ae ake 
SUBROUTINE REFRAC(JBEGIN, JEND,NPTS, IBEGIN, IEND, ISTART,M) 
COMMON/A/ C(60,20),RK(60,20),Y(60, 20) ,DEEP(60, 20) , ALPHAS(60, 20) 
COMMON/AA/YZERO(60) 
COMMON/B/ THETA(60,20),QXTOT(60), OLDANG(60,20), DY(60,20) 
COMMON/C/ H(60,20),CG(60, 20) ,HOLD(60, 20) ,HB(60, 20) , YB(60) 
COMMON/N USED/JUSE,T,CO,CGEN, CGGEN, ANGGEN, DX, BERM, THETAO( 10) , MMAX 
COMMON/D/SIGMA,G,ELO, JUMAX, IMAX,PI, TWOPI ,PI02,HGEN, IVET(10),SJUETTY 
COMMON/G/ IBREAK(60) , HNONBR( 20) 


COMMON/ZZZ/NTIME 

DIMENSION UJBEGIN(60), JEND(60) 
CRORE Ee Ae EEE ee THIS SUBROUTINE WILL DETERMINE H AND 
CR ek ESE ERE Ce THETA AT THE MID PT OF Y VALUES. 
C***TAU IS THE FACTOR WHICH RECOUPLES THE REFRACTION EQS.SEE ABBOTT 

TAU=0.25 


C*MUST PRESCRIBE THE WAVE ANGLE AT THE OUTERMOSTCONTOUR BOX 
C*SNELL’S LAW WILL BE USED TO START THINGS OFF. 
C*THETA(I,JU) WILL BE AT AREA’S CENTER AND WILL USE Y(I,JU) IN NEG Y-DIR 
C*WILL INITIALIZE ALL THETA’S USING SNELL‘S LAW. 
DO 206 I=IBEGIN, IEND 
C*INITIALIZE TWO J-VALUES BEYOND UJUMAX,IF IN REGION 1. 
IF (JEND(I).EQ.UMAX) JINIT=2 
IF (UEND(1I).NE.UJUMAX) JINIT=O 
DO 206 J=JBEGIN(I),JEND(I)+JINIT 
C*MUST CORRECT FOR THE CONTOUR ORIENTATION, ALPHAS. 
IF(I.NE.IBEGIN) GO TO 960 
ALPHAS(1,U)=ATAN((0.5*(Y(1+1,JU)+Y¥(14+1,U+1))-O.5*(Y(1, J) 
* +Y(I,U+1)))/DX) 
GO TO 962 


74 


65300 960 IF(I.NE.IEND) GO TO 961 


65400 ALPHAS(I,JU)=ATAN((0.5*(Y(1I,JU)+Y(1,JU+1))-O.5*(Y(1I-1,U) 
65500 * +Y(I-1,U+1)))/DX) 

65600 GO TO 962 

65700 961 ALPHAS(I,JU)=ATAN((0.5*(Y(1I+1,U)+Y¥(I+1,U+1))-0.5* 
65800 * (Y¥(I-1,U)+Y(I-1,U+1)))/(2.*DX)) 

65900 962 DALPHA=ANGGEN-ALPHAS(TI,uJ) 

66000 THETA(I,JU )=ARSIN((C(I,JU )/CGEN)*SIN(DALPHA) ) 

66100 C*MUST GET THETA WRT THE X-AXIS. 

66200 THETA(I,JU)=THETA(I,J)+ALPHAS(I,U) 

66300 206 CONTINUE 


66400 C*NOW, WE MUST COMP THE BOUN WAVE HTS SO THE HTS CAN BE COMPUTED. 
66500 CW HEE AUS Exit EE: Q)e sree As DEL DOT (E*CG)=0.0 


66600 C*NOW WE WILL CORRECT THE HT FOR SHOALING AND REFRACTION TO THE B.C. 
66700 C*WILL ALSO INITIALIZE H’S WITH THESE EQUATIONS FOR ENTIRE ARRAY. 
66800 DO 500 I=IBEGIN, IEND 

66900 C*INITIALIZE TWO JU-VALUES BEYOND UMAX IF IN REGION 1. 

67000 IF (JEND(I) .EQ.UMAX) JINIT=2 

67100 IF (JEND(I).NE.UMAX) JINIT=O 

67200 DO 500 J=JBEGIN(I1),JEND(I)+JUINIT 

67300 H(1I,JU) =HGEN*SQRT(CGGEN/CG(I,JU) )*SQRT(COS(ANGGEN) /COS(THETA(1, 
67400 * J))) 

67500 IF(HB(I,JU).LT.H(I,u)) H(1I,JU)=HB(I,uU) 

67600 500 CONTINUE 

67700 (OE oe ane ae eo oe ee ee a re Oe ST et aD le me iin ee 
67800 CII I I IIR RO III III IOI III I FR FO 
67900 C*LET’S FILL THE DY ARRAY. 

68000 C*DY WILL BE INDEXED AS THE THETA TO WHICH WE ARE GOING. 

68100 DO 209 I=IBEGIN, IEND 

68200 DO 209 J=JBEGIN(I )+1,JEND(I) 

68300 DY(I,JU-1)=0.5*(Y(1I,JU-1)+Y(1I,J))-O.5*(Y(I,JU)+Y¥(1,U+1)) 

68400 209 CONTINUE 

68500 NITERS=100 

68600 DO 100 NITER=1,NITERS 

68700 SUMANG=0.O 

68800 C*DO "6O LOOP" GOES FROM 2 TO IMAX IF ISTART =IBEGIN 

68900 c*DO "60 LOOP" GOES FROM IMAX-1 TO 1 IF ISTART=IEND 

69000 DO 6O II=IBEGIN, IEND 

69100 C*MUST HAVE IT SET UP SO THAT THE KNOWN BOUNDARIES ANGLES AREN’T RECOMP 
69200 IF(ISTART.EQ. IBEGIN) I=I1 

69300 IF(ISTART.EQ.IBEGIN .AND. I1.E£Q.IBEGIN) GO TO 60 

69400 IF(ISTART.EQ.IEND) I=IEND-II+IBEGIN 

69500 IF(ISTART.EQ.IEND .AND. I1.EQ.IEND) GO TO 60 

69600 C*ADX EQUALS ACTUAL DELTA X ACROSS SPACE STEP. 

69700 C*ONLY ON BOUNDARIES WHERE FORWARD OR BACKWARD DIFFERENCING. 

69800 IF(I.NE.IBEGIN) GO TO 6 

69900 ADX=DX 

70000 IP=I+1 

70100 IM=I 

70200 GO TO 12 

70300 6 IF(I.NE.IEND) GO TO 10 

70400 ADX=DX 

70500 IP=I 

70600 IM=I-1 

70700 GO TO 12 

70800 10 ADX=2 .O*DX 

70900 IP=I+1 

71000 IM=I-1 

71100 12 CONTINUE 

71200 DO 40 J=JBEGIN(I),JEND(I)-1 

71300 C*WILL GO FROM (JMAX-1) TO 1 BECAUSE THAT’S THE DIR WAVE COMES IN FROM. 
71400 JJ=JEND(I )-1-JU+JUBEGIN(T) 

71500 OLDANG(I,JUJU)=THETA(I, Ju) 

71600 C*LOCATE MIDPOINT BETWEEN TWO ADJACENT BLOCK CENTERS 


71700 C*BECAUSE THETA’S JJ-VALUE IS THE SAME AS THE FIRST SHOREWARD Y VALUE 
71800 C*MUST. USE JJ, UJ+1, AND JJ+2 TO COMPUTE YBAR. 


71900 YBAR=0. 25*(Y(I,Ju)+2.0*Y(1I,JJU+1)+Y(1, UU+2) ) 
72000 C*LOCATE APPROPRIATE INDICES ON IP AND IM GRID LINES. 
72100 IMINUS=- 1 

72200 IPLUS=+1 

72300 CALL LOC(IM,JJ,JOIM, JSIM, YBAR, IMINUS) 

72400 CALL LOC(IP,UJ,JOIP,JSIP,YBAR,IPLUS) 


72500 C*NOW USE THE CONSERVATION OF WAVES EQUATION.............. 
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72600 
72700 
72800 
72900 
73000 
73100 
73200 
73300 
73400 
73500 
73600 
73700 
73800 
73900 
74000 
74100 
74200 
74300 
74400 
74500 
74600 
74700 
74800 
74900 
75000 
75100 
75200 
75300 
75400 
75500 
75600 
75700 
75800 
75900 
76000 
76100 
76200 
76300 
76400 
76500 
76600 
76700 
76800 
76900 
77000 
77100 
77200 
77300 
77400 
77500 
77600 
77700 
77800 
77300 
78000 
78100 
78200 
78300 
78400 
78500 
78600 
78700 
78800 
78900 
79000 
79190 
79200 
79300 
79400 
79500 
79600 
79700 
79800 


PART 1C=RK(I, JU+1)*SIN(THETA(I, JU+1)) 
PART2=-DY(I,JUJ)/ADX 
C*WILL LINEARLY INTERPOLATE TO DETERMINE RK*COS(THETA) AT I+1 AND I-1. 
C*IF NO ADJ SHOREWARD PT EXISTS, PUT IN ZERO FOR TERMS IN GOV. EQ. 
IF(USIM.NE.O) GO TO 301 
PART3B=0.0 
GO TO 302 
301 TOPIM=RK( IM, JOIM- 1)*COS(THETA(IM, JOIM-1)) 
BOTIM=RK( IM, JSIM) *COS(THETA(1IM,USIM) ) 
TOTALB=0.5*(Y( IM, JOIM)+Y( IM, JOIM-1))-0.5*(Y( IM, USIM+1)+Y(IM,JSIM) ) 
DUMB=0.5*(Y( IM, JOIM)+Y(IM, JOIM-1))-YBAR 
PART3B=((TOTALB-DUMB ) *(TOPIM-BOTIM)/TOTALB)+BOTIM 
302 IF(USIP.NE.O) GO TO 303 
PART3A=0.0 
GO TO 304 
303 TOPIP=RK(IP,UOIP-1)*COS(THETA( IP, JOIP-1)) 
BOTIP=RK(IP,JUSIP)*COS(THETA(IP,USIP) ) 
TOTALA=0.5*(Y(IP,JOIP)+Y(IP,JOIP-1))-0.5*(Y(IP,USIP+1)+Y(IP,USIP) ) 
DUMA=0.5*(Y(IP,JOIP)+Y(IP,JOIP-1))-YBAR 
PART3A=(( TOTALA-DUMA ) *(TOPIP-BOTIP)/TOTALA)+BOTIP 
304 PART3=PART3A-PART3B 
C*NOW MUST FIND RK*SIN(THETA) FOR I+1 AND I-14 AT J+1 
YBARP=O. 25*(Y(1I,JJU+1)+2.*Y(1,JJU+2)+Y(I,JU+3)) 
CALL LOC(IM, JJ+1,JPOIM, JUPSIM, YBARP, IMINUS) 
CALL LOC(IP,JJU+1,JUPOIP,UPSIP, YBARP,IPLUS) 
IF(UPSIM.NE.O) GO TO 305 
PART 1B=0.0 
GO TO 306 
305 TOPM=RK( IM, JPOIM-1)*SIN( THETA( IM, JUPOIM-1)) 
BOTM=RK(IM, JPSIM)*SIN( THETA( IM, JUPSIM) ) 
TOTB=0.5*(Y( IM, JPOIM)+Y(IM, JPOIM-1))-0.5*(Y( IM, JPSIM+1)+ 
* Y( IM, JPSIM) ) 
DUMPB=0.5*(Y (IM, JPOIM)+Y(IM, JPOIM-1))-YBARP 
PART 1B=(( TOTB-DUMPB ) * (TOPM-BOTM)/TOTB)+BOTM 
306 IF(UPSIP.NE.O) GO TO 307 
PART 1A=0.0 
GO TO 308 
307 TOPP=RK(IP,JUPOIP-1)*SIN(THETA(IP,JPOIP-1)) 
BOTP=RK(IP,JUPSIP)*SIN(THETA(IP,UPSIP) ) 
TOTA=0.5*(Y(IP,JUPOIP)+Y(IP,JUPOIP-1))-0.5*(Y(IP,UPSIP+1)+Y(IP,JUPSIP 
ce )) 
DUMPA=0.5*(Y(IP,UPOIP)+Y(IP,JPOIP-1))-YBARP 
PART 1A=(( TOTA-DUMPA ) *( TOPP-BOTP )/TOTA)+BOTP 
308 PART 1=TAU*PART1B+(1.-2.*TAU)*PART1C+TAU*PART1A 
IF(UPSIM.EQ.0O)PART1=(1.-TAU) *PART1C+TAU*PARTI1A 
IF(UPSIP.EQ.0O)PART 1=TAU*PART 1B+(1.-TAU) *PART1C 
ARG=( (PART 1+PART2*PART3)/RK(I, Ju) ) 
C*IF THE ROUTINE IS TO BLOWUP,USE SNELLS LAW. 
IF(ABS(ARG).LE.1.0) GO TO 41 
ARG=(C(1I,JUJU)/C(1I, Ju+1))*SIN(THETA(I,JuU+1)) 
IF(ARG.GT.1.0) ARG=1.0 
THETA(1I,JUJU)=ARSIN(ARG) 
GO TO 42 
41 THETA(I,JJ)=ARSIN( ARG) 
42 THETA(I,JJU)=0.5*(THETA(I,JJ)+OLDANG(I, Ju) ) 
SUMANG=SUMANG+(ABS(THETA(I, JJ) -OLDANG(I,uUU))) 
40 CONTINUE 
60 CONT INUE . 
C*MUST EJECT IF WE HAVE REACHED AN ACCEPTABLE ITERATION ERROR 
C*IF THE SUM OF THE ABSOLUTE VALUE OF ANGLE CHANGES DURING AN ITERATION 
c* AVERAGES LESS THAN 0.02 DEGREES PER GRID ITS CLOSE ENOUGH. 
IF (SUMANG.LT. (NPTS*O.0035) ) GO TO 215 
IF(NITER.GE.50) GO TO 215 
100 CONTINUE 
WRITE(6, 803) 
215 CONTINUE 
C*ITERATION LOOP FOR THE WAVE HEIGHT. 
DO 501 NITER=1,NITERS 
SUMH=0.0 
DO 510 II=IBEGIN,IEND 
C*MUST HAVE IT SET UP SO THAT THE KNOWN BOUNDARIES HTS. AREN’T RECOMP 
IF(ISTART.EQ.IBEGIN) I=II 
IF(ISTART.EQ.IBEGIN .AND. I1.EQ.IBEGIN) GO TO 510 
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79900 
80000 
80100 
80200 
80300 
80400 
80500 
80600 
80700 
80800 
80900 
81000 
81100 
81200 
81300 
81400 
81500 
81600 
81700 
81800 
81900 
82000 
82100 
82200 
82300 
82400 
82500 
82600 
82700 
82800 
82900 
83000 
83100 
83200 
83300 
83400 
83500 
83600 
83700 
83800 
83900 
84000 
84100 
84200 
84300 
84400 
84500 
84600 
84700 
84800 
84900 
85000 
85100 
85200 
85300 
85400 
85500 
85600 
85700 
85800 
85900 
86000 
86100 
86200 
86300 
86400 
86500 
86600 
86700 
86800 
86900 
87000 


IF (ISTART.EQ. TEND) I=IEND-II+IBEGIN 
IF(ISTART.EQ.IEND .AND. I.EQ.IEND) GO TO 510 


C*ADX EQUALS ACTUAL DELTA X ACROSS SPACE STEP. 
C*ONLY ON BOUNDARIES WHERE FORWARD OR BACKWARD DIFFERENCING. 


503 


504 


505 


312 


313 


315 


316 


317 


318 


IF(I.NE.IBEGIN) GO TO 503 
ADX=DX 
IP=I+1 
IM=I 
GO TO 505 

IF(I.NE.IEND) GO TO 504 
ADX=DX 
IP=I 
IM=I-14 
GO TO 505 

ADX=2.0*DX 
IP=I+1 
IM=I-1 

CONT INUE 
DO 502 J=JBEGIN(I),JEND(I)-1 
JJ=JEND(I)- 1-JU+UBEGIN(I) 
HOLD(I,JUJU)=H(I,uu) 
YBAR=O.25*(Y(1I,JU)+2.0*Y(1, JU+1)+Y(1,JJ+2) ) 
CALL LOC(IM, JJ,JOIM,JUSIM, YBAR, IMINUS) 
CALL LOC(IP,UJU,JOIP,USIP, YBAR, IPLUS) 
PART13=(H(I,JJU+1)**2.)*CG(I,JJ+1)*COS(THETA(I, JuU+1)) 
PART2=DY(I,JUJU)/ADX 
IF(JUSIM.NE.O) GO TO 311 
PART4B=0.0 
GO TO 312 
TOPIMH=(H( IM, JOIM-1)**2. )*CG( IM, JOIM-1)*(SIN(THETA( IM, JOIM-1))) 
BOTIMH=(H( IM, JSIM) **2. )*CG( IM, JSIM) *SIN( THETA( IM, JSIM) ) 
TOTALB=0.5*(Y( IM, JOIM)+Y (IM, JOIM-1))-0.5*(Y( IM, JSIM+1)+Y(IM,JUSIM) ) 
DUMB=0.5*(Y( IM, JOIM)+Y(IM, JOIM-1))-YBAR 
PART4B=((TOTALB-DUMB ) * (TOP IMH-BOTIMH)/TOTALB )+BOTIMH 
IF(JUSIP.NE.O) GO TO 313 
PART4A=0.0 
GO TO 314 
TOPIPH=(H(IP,JOIP-1)**2.)*CG(IP, JOIP-1)*SIN(THETA(IP,JOIP-1)) 
BOTIPH=(H(IP, JSIP)}**2.)*CG(IP,JSIP)*SIN(THETA(IP,USIP)) 
TOTALA=0.5*(Y(IP,JOIP)+Y(IP,JOIP-1))-0.5*(Y(IP,JUSIP+1)+Y(IP,USIP) ) 
DUMA=0.5*(Y(IP,JOIP)+Y(IP,JOIP-1))-YBAR 
PART 4A=((TOTALA-DUMA ) *(TOPIPH-BOTIPH)/TOTALA )+BOTIPH 
PART4=PART4A-PART4B 
YBARP=O0.25*(Y(I,JU+1)+2.*Y(1,JUJU+2)+Y(1, JUU+3)) 
CALL LOC(IM, JJ+1,JPOIM, UPSIM, YBARP, IMINUS) 
CALL LOC(IP,JJ+1,JPOIP,UPSIP,YBARP, IPLUS) 
IF(UPSIM.NE.O) GO TO 315 
PART12=0.0 
GO TO 316 
TOPMH=(H( IM, JPOIM-1)**2)*CG( IM, JPOIM-1)*COS(THETA(IM, JPOIM-1) ) 
BOTMH=(H( IM, UPSIM) **2)*CG( IM, UPSIM) *COS(THETA( IM, JUPSIM) ) 
TOTB=.5*(Y( IM, JPOIM)+Y(IM, JPOIM-1))-.5*(Y( IM, JPSIM+1)+Y(IM,JPSIM) ) 
DUMPB=0.5*(Y(IM, JPOIM)+Y (IM, JPOIM-1))-YBARP 
PART 12=((TOTB-DUMPB ) * (TOPMH-BOTMH) /TOTB )+BOTMH 
IF(UPSIP.NE.O) GO TO 317 
PART11=0.0 
GO TO 318 
TOPPH=(H( IP, UPOIP-1)**2)*CG(IP, JUPOIP-1)*COS(THETA(IP,UPOIP-1)) 
BOTPH=(H( IP, UPSIP)**2)*CG(IP,UPSIP)*COS(THETA(IP,UPSIP) ) 
TOTA=.5*(Y(IP, UPOIP)+Y(IP, UPOIP-1))-.5*(Y(IP,JUPSIP+1)+Y(IP,UPSIP) ) 
DUMPA=0.5*(Y(IP,JUPOIP)+Y(IP,JUPOIP-1))-YBARP 
PART 11=((TOTA-DUMPA ) *(TOPPH-BOTPH) /TOTA)+BOTPH 
PART 1H=TAU*PART 12+(1.-2.*TAU) *PART13+TAU*PART 11 
IF(UPSIM.EQ.O)PART1H=(1.-TAU) *PART13+TAU*PART 11 
IF(JUPSIP.EQ.0)PART 1H=TAU*PART 12+(1.-TAU)*PART13 
ARG=( (PART 1H+PART2*PART4)/(CG(I,JUJU)*COS(THETA(I,uJ)))) 


C*IF THERE IS TO BE AN INVALID SQRT,USE LINEAR SHOALING. 


IF(ARG.GE.O.) GO TO 44 

ARG=(CG(I, JU+1)*COS(THETA(I,JU+1)))/(CG(1, JU) *COS(THETA(I, JJ) )) 
IF(ARG.LT.0O.0) ARG=0.0 

H( 1, UJ) =H(1I,JJ+1)*SQRT (ARG) 

GO TO 45 
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H(1I,JU)=SQRT (ARG) 
H(1I,JJU)=0.5*(H(1,JU)+HOLD(I, UU) ) 
HNONBR( JJ) =H(I, JJ) 


C*IBREAK(I)=JJU, THEREFORE JJ WILL BE LEEWARD SIDE OF GRID AT INIT BREAK 


501 
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* 


IF(HB(I,JU) .LT. H(1,JJ) .AND. HB(1I,JJU+1).GE.HNONBR(JuU+1) ) 
IBREAK(I)=Ju 

IF(HB(I,JJU).LT.H(I,uu)) H(1,JUJU)=HB(1, JJ) 

SUMH=SUMH+ABS(H(I,JJU)-HOLD(I, Ju) ) 

CONT INUE 

CONTINUE 

IBREAK( IEND)=IBREAK(IEND-1) 

IBREAK(IBEGIN) =IBREAK( IBEGIN+1) 

IF(SUMH.LT. (NPTS*O.01)) GO TO 507 

IF(NITER.GE.50) GO TO 507 

CONTINUE 

WRITE(6,803) 

CONT INUE 

FORMAT(2X,4(F15.5),////) 

FORMAT(2X, "AFTER NITERS ITERATIONS, CONVERGENCE WAS NOT REACHED") 


FORMAT(2X,"THE WAVE HT. ROUTINE CONVERGED IN, NITER= ",15,//) 
FORMAT(2X, "THIS IS MY CHECKING WRITE STATEMENT") 

FORMAT(2X,"THE WAVE ANGLE ROUTINE CONVERGED IN, NITER= ",15,//) 
RETURN 

END 


CI RI OR II I RO I III ICI III I II II II ICI II I IOI I a tO Ok oe 


SUBROUTINE DIFF(RHOND, THETAO, ANGLE, AMP) 


C****DIFFRACTION ABOUT SEMI INFINITE BREAKWATER (PENNEY-PRICE ) 


PI=3.14159265 
ABSS=SIN(O.5*(ANGLE-THETAO) ) 
ABSP=SIN(O.5*(ANGLE+THETAO) ) 
ABC=COS(ANGLE-THETAOQ) 

ABC 1=COS(ANGLE+THETAO) 

XX =RHOND* ABC 

XXC=COS(XX) 

XXS=SIN(XX) 

XX 1=RHOND*ABC 1 

XXC1=COS(XX1) 

XXS1=SIN(XX1) 
AL=SQRT(RHOND/PI ) 
SIG=2.0*AL*ABSS 
SIGP=-2.0*AL*ABSP 

CALL FRES(SIG,C,S,FR,FI) 

CALL FRES(SIGP,CP,SP,FRP,FIP) 
SUM1=XXC*FR+XXS*FI+XXC1*FRP+XXS1*F IP 
SUM2=XXC*FI-XXS*FR+XXC1*FIP-XXS1*FRP 
AMP=SQRT (SUM1**2+SUM2**2) 
RETURN 

END 
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SUBROUTINE FRES(A,C,S,FR,FI) 


C*FRESNEL INTEGRAL SUBROUTINE****AFTER ABROMOWITZ AND STEGUN. 
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Z=ABS(A) 

PO2=1.5707963 
FZ=(1.0+0.926*Z)/(2.0+1.792*Z+3. 104*Z*Z) 
GZ=1.0/(2.0+4. 142*2+3.492*Z*Z2+6 .670*Z*Z*Z) 
XX =PO2*Z*Z 

CZ=COS(XX) 

SZ=SIN(XX) 

C=0.5-GZ*CZ+FZ*SZ 

S=O0.5-FZ*CZ-GZ*SZ 

IF(A.GT.O.0) GO TO 50 

C=-C 

S=-S 

FR=0.5*(1.0+C+S) 

FI=-0.5*(S-C) 

RETURN 

END 
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SUBROUTINE PREDIF 

COMMON/A/ C(60,20),RK(60,20),Y(60, 20) ,DEEP(60, 20) , ALPHAS(60, 20) 
COMMON/AA/YZERO(60) 

COMMON/B/ THETA(60,20),QXTOT(60), OLDANG(60,20), DY(60,20) 
COMMON/C/ H(60,20),CG(60, 20) ,HOLD(60, 20) ,HB(60, 20), YB(60) 
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COMMON/N USED/JUSE,T,CO,CGEN,CGGEN, ANGGEN, DX’, BERM, THETAO( 10) ,MMAX 
COMMON/D/SIGMA,G,ELO,JUMAX, IMAX,PI, TWOPI ,PI02,HGEN, IJET(10),SUETTY 
COMMON/G/IBREAK(60) , HNONBR( 20) 
DIMENSION J1(60),U2(60),J1REF(60) ,J3REF(60) 
C*THIS SUB CALCS WHERE DIFFRACTION GOVERNS AND WHERE REFRACT GOVERNS. 
C*IT WILL CALL REFRAC FOR OFFSHORE AREA(OFF TIP OF STRUCTURE). 
C*THEN IT WILL DO THE SHADOW ZONE USING DIFF(IF THETAO .NE.O.O) 
C* IT WILL THEN FINISH THE OTHERS USING REFRAC AGAIN. 
C*LET’S ZERO-OUT THE DIMENSIONED ARRAYS. 
DO 1000 I=1,IMAX+2 


Ji(1)=0.0 
J2(1)=0.0 
J1REF(1I)=0.0 
1000 JU3REF(1I)=0.0 
C*NOW, LETS FIND C,CG,RK,HB, AND WVNUM. 


DO 202 I=1,IMAX 
DO 202 J=1,JUMAX+2 
DEPTH=DEEP(I,u) 
CALL WVNUM(DEPTH, T,DUMK) 
RK( I,J) =DUMK 
C(1I,JU)=CO*TANH(RK(1I,JU)*DEEP(I,U)) 
EN=0.5*(71.0+((2.*RK(1I,JU)*DEEP(I,U))/SINH(2.*RK(I,JU)*DEEP(I,J)))) 
CG(I,JU)=EN*C(I,U) 
HB(I,J)=0.78*DEEP(I,uJ) 
202 CONTINUE 
C*WILL ATTRIB AN EQUAL REACH TO EACH SIDE OF EACH M-GROIN. 
DO 200 M=1,MMAX 
IDUML=1 
IF(M.NE.1) IDUML=(IJET(M)+IJET(M-1))/2 
IDUMR=IMAX 
IF(M.NE.MMAX) IDUMR=(IJET(M)+IJET(M+1))/2 
NPTS=O 
DO 1 I=IDUML,IDUMR 
DO 2 J=1,JMAX 
IF(Y(1I,J).LT.SJETTY) GO TO 14 
Ji(1I)=J 
J2(1)=JMAX 
GO TO 15 
144 CONTINUE 
2 CONTINUE 
145 CONTINUE 
C*IF NO STRUCT IS PRESENT(SJETTY=0.0), DO REFRAC THRUOUT GRID SYSTEM 
IF(SJETTY.EQ.0.0) J1i(I)=1 
4 CONTINUE 
DO 16 I=IDUML,IDUMR 
C* ‘REFRAC’ STARTS ON THE NEXT TO LAST JU-CONTOUR,NOT THE LAST! 
DO 16 J=J1(1),J2(1)-1 
16 NPTS=NPTS+14 
C*WILL NOW DO THE REFRACT FOR’THE REGION 1 AREA. 
C*ISTART REPRESENTS THE DIRECTION THE SWEEPS WILL BEGIN FROM. 
C*WILL USE DUMMY IMAX,IJUET,IJET+1 IN CALL STTS SO IBEGIN,IEND, AND 
C***ISTART WON’T CHANGE THEM.MUST RESET AFTER EACH CALL REFRAC. 
IMAXT=IDUMR 
IJETT=IJET(M) 
IJETP1=IJET(M)+1 
IDUMLL=IDUML 
IF(ANGGEN.GE.O.0) CALL REFRAC(JU1,J2,NPTS, IDUMLL, IMAXT, IDUMLL,M) 
IF(ANGGEN.LT.O.0) CALL REFRAC(Ji1,J2,NPTS, IDUMLL, IMAXT, IMAXT,M) 
IMAXT=IDUMR 
IJUETT=IJET(M) 
IJETP1=IJET(M)+1 
IDUMLL=IDUML 
JDUMN=J4(IJET(M) ) 
JDUMS=JU1(IJET(M)+1) 
XDISTN=(IJUET(M)-1.0)*DX+DX/2. 
ELTIP=T*0.5*(C(IJET(M) , JDUMN)+C(IJUET(M)+1,JDUMS) ) 
C*NOW MUST CHECK THE ANGLE AT THE STRUCTURE’S TIP TO SEE WHERE SHAD ZONE 
C*IF NO STRUCT PRESENT(SJETTY=0.0), FUTHER REFRAC/DIFF UNNECESSARY. 
IF(SUJETTY.EQ.0.0) GO TO 13 
THETAO(M)=0.5*(THETA(IJET(M), JDUMN)+THETA( IJVET(M)+1,JUDUMS) ) 
HINC=0.5*(H(IJET(M) , JDUMN)+H( IJET(M)+1,JUDUMS) ) 
IF(THETAO(M))10,11,12 
C*THIS SECTION HANDLES REFRAC/DIFF IF THETAO<O.O. 


79 


101700 10 CONTINUE 


101800 C*FIRST ALL OF REGION 2 WILL GET REFRACTED. 

101900 NPTS=O 

102000 DO 100 I=IJET(M)+1,IDUMR 

102100 J2(1)=JU1(1) 

102200 100 Ji(1I)=1 

102300 DO 101 I=IJET(M)+1,IDUMR 

102400 DO 101 J=JU1(1),J2(I)-1 

102500 101 NPTS=NPTS+1 

102600 IMAXT=IDUMR 

102700 IDUMLL=IDUML 

402800 IJETT=IJET(M) 

1402900 IJETP1=IJET(M) +1 

103000 CALL REFRAC(JU1,JU2,NPTS,IJETP1, IMAXT, IMAXT,M) 
103100 IMAXT=IDUMR 

103200 IJETT=IJET(M) 

103300 IJETP1=IJET(M)+1 

103400 IDUMLL=IDUML 

103500 C*NOW MUST DO REGION 3 OF NEG THETAO CASE-SHADOW ZONE. 
103600 DO 102 I=IDUML,IJET(M) 

103700 J2(1)=J1(1) 

103800 102 J1i(1)=1 

103900 DO 103 I=IDUML,IJET(M) 

104000 JIREF(1)=1 

104100 DO 104 J=J1(1I),J2(1)+1 

104200 XCOOR=(I-1.0)*DX 

104300 YCOOR=0.5*(Y(1I,U)+Y¥(1,U+1)) 

104400 ANGLE=ATAN( (XDISTN-XCOOR)/(SJETTY-YCOOR) ) 
104500 IF(YCOOR.GT.SJETTY) ANGLE=PI+ANGLE 

104600 C*IF MOST SHOREWARD PT OUT OF SHAD ZONE, SO ARE THE OTHERS FOR THAT I. 
104700 IF (ABS(ANGLE).GT.ABS(THETAO(M) ) ) GO TO 105 
104800 RAD=SQRT((XDISTN-XCOOR) **2+(SJETTY-YCOOR) **2) 
104900 RHOND=RAD*TWOPI/ELTIP 

105000 C*DIFFRACTION TREATS THE POS THETAO CASE. 

105100 THE =ABS(THETAO(M) ) 

105200 CALL DIFF (RHOND, THE,ANGLE, AMP) 

105300 H( I,J) =AMP*HINC 

105400 ANGRAD=-ANGLE 


105500 C*WILL NOW REFRACT DIFF WAVES IN THE SHAD ZONE USING SNELL’S. 
105600 CTIP=ELTIP/T 


105700 ALPHAS(I1,J)=ATAN((0.5*(Y(I+1,JU)+Y¥(1I+4,JU+1))-O0.5* 
105800 * (Y(I-1,JU)+Y(1I-1,JU+1)))/(2. *DXx)) 

105300 IF(I.EQ.IJET(M) )ALPHAS(1,JU)=ATAN((0.5*(Y(1I,JU)+Y(1I,U+1))-0.5*(Y(I-1 
106000 * ,U)+¥(1-1,JU+1)))/DX) 

1406100 DALPHA=ANGRAD-ALPHAS(I,U) 

106200 THETA(1I,JU)=ARSIN((C(1,JU)/CTIP)*SIN(DALPHA) ) 
106300 THETA(1I,JU)=THETA(I,JU)+ALPHAS(I,uJ) 

106400 C*MUST CHECK TO SEE IF WAVE WOULD HAVE BROKEN. 

106500 IF(HB(I,J).LE.H(1I,JU).AND.HB(1,JU+1).GT.H(1,U+1))IBREAK(I)=J 
106600 IF(HB(I,JU).LT.H(I,u)) H(1,JU)=HB(I,U) 

106700 104 CONTINUE 

106800 GO TO 103 

106900 105 J1REF(1I)=uJ 

107000 103 CONTINUE 

107100 C*NOW MUST DO REFRACTION FOR REGION 4. 

107200 NPTS=O 

1407300 DO 106 I=IDUML,IJET(M) 

107400 DO 106 U=J1REF(I),J2(1I)-1 

1407500 106 NPTS=NPTS+1 

107600 IDUMLL=IDUML 

107700 IMAXT=IDUMR 

107800 IJETT=IJET(M) 

107900 IJETP1=IJET(M)+1 

108000 CALL REFRAC(J1REF,J2,NPTS, IDUMLL, IJETT, IDUMLL,M) 
108100 IDUMLL=IDUML 

108200 IMAXT=IDUMR 

108300 IJETT=IVJET(M) 

408400 IJETP1=IJET(M)+14 

108500 GO TO 13 

408600 C*THIS HANDLES REFRAC/DIFF IF THETAO IS O.O. 

108700 C*FOR THIS CASE, ONLY THREE REGIONS EXIST. 

108800 44 CONTINUE 

108900 NPTS=O 
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109000 DO 120 I=IDUML,IVJET(M) 


109100 J2(1)=JU1(1) 

109200 120 Ji(I)=1 

109300 DO 121 I=IDUML,IVET(M) 
109400 DO 1241 J=J1(1I),J2(1I)-1 
109500 421 NPTS=NPTS+1 

109600 IMAXT=IDUMR 

108700 IDUMLL=IDUML 

109800 IJETT=IJET(M) 

109900 IJETP1=IJET(M) +4 

110000 CALL REFRAC(Ji,J2,NPTS, IDUMLL, IVETT, IDUMLL,M) 
110100 IMAXT=IDUMR 

110200 IJETT=IJET(M) 

110300 IJETP1=IJET(M) +1 

110400 IDUMLL=IDUML 

110500 DO 122 I=IJVET(M)+1,IDUMR 
110600 J2(1I)=U1(1) 

410700 122 Ji(I)=1 

110800 NPTS=O 

110900 DO 123 I=IJET(M)+1, IDUMR 
111000 DO 123 J=Ji(I),J2(1)-1 
111100 1423 NPTS=NPTS+1 

111200 IMAXT=IDUMR 

111300 IDUMLL=IDUML 

111400 IJETT=IVJET(M) 

111500 IJETP1=IJET(M)+1 

111600 CALL REFRAC(JU1,JU2,NPTS,IJVETP1, IMAXT, IMAXT,M) 
111700 IMAXT=IDUMR 

111800 IJETT=IJET(M) 

111900 IJETP1=IJET(M)+4 

112000 IDUMLL=IDUML 

112100 GO TO 13 


112200 C*THIS SECTION HANDLES REFRACT/DIFF IF THETAO>0.0 
112300 12 CONTINUE 
112400 C*FIRST, REGION 2- ALL REFRACTION. 


112500 NPTS=O 

112600 DO 110 I=IDUML,IJET(M) 

112700 J2(1)=J1(1) 

112800 110 Ji(1)=1 

112900 DO 111 I=IDUML,IJET(M) 

113000 DO 111 J=J1(I),J2(1)-1 

113100 111 NPTS=NPTS+1 

113200 IMAXT=IDUMR 

113300 IDUMLL=IDUML 

113400 IJETT=IJET(M) 

113500 IJETP1=IJET(M)+1 

113600 CALL REFRAC(JU1,JU2,NPTS, IDUMLL, IJETT, IDUMLL,M) 
113700 IMAXT=IDUMR 

113800 IJETT=IJET(M) 

113900 IJETP1=IJET(M)+14 

114000 IDUMLL =IDUML 

114100 C*NOW WILL DO REGION 3 OF THE POS THETAO CASE. 
114200 DO 112 I=IJET(M)+i,IDUMR 

114300 J2(1I)=Ji(1) 

114400 112 vi(1)=1 

114500 DO 113 I=IJET(M)+1,IDUMR 

1414600 JIREF(1)=1 

114700 C*WILL GO ONE PT. BEYOND JU2(I) TO MAKE SURE OUTOF DIFF ZONE. 
114800 DO 114 J=Ji(I),J2(1)+14 

114900 XCOOR=(I-1.0)*DX 

115000 YCOOR=0.5*(Y(1,JU)+Y(1,uU+1)) 

115100 ANGLE =ATAN( (XCOOR-XDISTN)/(SJETTY-YCOOR) ) 
115200 IF(YCOOR.GT.SUJETTY) ANGLE =PI+ANGLE 

115300 C*IF LEAST J-VALUE IS OUT OF SHAD ZONE,SO ARE OTHER JU’S.(FOR EACH I) 
415400 IF (ANGLE .GT.ABS(THETAO(M) ) ) GO TO 115 
115500 RAD=SQRT( ( XCOOR-XDISTN) **2+(SJETTY-YCOOR) **2) 
415600 RHOND=RAD*TWOPI/ELTIP 

415700 THE=THETAO(M) 

115800 CALL DIFF(RHOND, THE,ANGLE, AMP) 

115900 ANGRAD=ANGLE 


116000 C*WILL NOW REFRACT DIFFRACTED WAVES IN SHAD ZONE USING SNELL"S. 
116100 CTIP=ELTIP/T 


116200 ALPHAS(1,JU)=ATAN((0.5*(Y(I+1,JU)+Y¥(1I+1,JU+1))-0.5* 
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= AGES o RAK =). 5 eea)) ) )/ C2.) ) 
IF(I.EQ. IUJET(M)+1)ALPHAS(I,JU)=ATAN((0.5*(Y(I+1,JU)+Y(1+1,U+1))-O.5* 
* (Y(I,JU)+Y¥(1I,J+1)))/DX) 
DALPHA=ANGRAD-ALPHAS(I,U) 
THETA(I,JU)=ARSIN((C(I,JU)/CTIP)*SIN(DALPHA) ) 
THETA(I,JU)=THETA(I,JU)+ALPHAS(I,uU) 
H(1,J)=HINC*AMP 
C*MUST CHECK TO SEE IF WAVE WOULD HAVE BROKEN. 
IF(HB(I,JU).LE.H(1,JU).AND.HB(1I,JU+1).GT.H(1I,JU+1))IBREAK(I)=u 
IF(HB(I,JU).LT.H(I,U)) H(I,J)=HB(I,J) 
114 CONTINUE 
GO TO 113 
115 JIREF(I)=J 
413 CONTINUE 
C*NOW MUST DO REFRAC FOR REGION 4. 
NPTS=O 
DO 116 I=IJET(M)+1,IDUMR 
DO 116 J=J1REF(I),J2(1)-1 
116 NPTS=NPTS+14 
IMAXT=IDUMR 
IDUMLL=IDUML 
IJETT=IJET(M) 
IJETP1=IJET(M)+1 
CALL REFRAC(J1REF,JU2,NPTS, IJETP1, IMAXT, IMAXT,M) 
IMAXT=IDUMR 
IJETT=IJET(M) 
IJVETP1=IJET(M)+1 
IDUMLL=IDUML 
143 CONTINUE 
200 CONTINUE 
RETURN 
END 
CE ee te OK ee ke fe ek 2 ee ke oe ok ok kok a ok kK ok ok oe ke ok kK ek kok 2K kk ok ek 2k ok ok oe ok 
SUBROUTINE LOC(IM,JJ,JOIM,JSIM, YBAR, IDUM) 
COMMON/A/ C(60,20),RK(60,20),Y(60, 20) ,DEEP(60, 20) , ALPHAS(60, 20) 
COMMON/AA/YZERO(60) 
COMMON/B/ THETA(60,20),QXTOT(60), OLDANG(60,20), DY(60,20) 
COMMON/C/ H(60,20),CG(60, 20) ,HOLD(60, 20) ,HB(60, 20), YB(60) 
COMMON/N USED/JUSE,T,CO,CGEN, CGGEN, ANGGEN, DX ,BERM, THETAO( 10) , MMAX 
COMMON/D/SIGMA,G,ELO, JMAX, IMAX,PI, TWOPI,PI02,HGEN, IJET(10),SUETTY 
C*SUBROUTINE LOC FINDS J-VALUES WHICH ARE GREATER AND LESS THAN YBAR. 
JOIM=2 
2 AA=0.5*(Y( IM, JOIM)+Y (IM, JOIM-1) ) 
IF(AA.GT.YBAR) GO TO 4 
JOIM=JOIM+ 1 
C*THE FOLLOWING IS REQ’D SO THAT DY/DX>0.5 
C*WILL DTERMINE K SIN THETA ON IM-LINE AT A DIST YBAR. 
C*WILL CALL THIS POINT JUSE+1 
IF(JOIM.LE.JUSE) GO TO 2 
JOIM=JUSE+1 
Y( IM, JOIM)=YBAR 
C* DEPTH AT THIS POINT WILL BE COMP ASSUMING CONST BEACH SLOPE ON I=IM 
DEL=.5*(Y( IM, JOIM-1)+Y( IM, JOIM-2))-.5*(Y(IM, JOIM-2)+Y(IM,JOIM-3) ) 
BSLOPE=(DEEP(IM, JOIM-2)-DEEP(IM, JOIM-3))/DEL 
DEEP( IM, JOIM-1)=DEEP(IM, JOIM-2)+BSLOPE*(Y(IM, JOIM)-Y(IM, JOIM-1)) 
DEPTH=DEEP(IM, JOIM- 1) 
CALL WVNUM(DEPTH,T,DUMK) 
RK( IM, JOIM- 1) =DUMK 
C( IM, JOIM-1)=CO*TANH(RK( IM, JOIM-1)*DEEP( IM, JOIM-1)) 
EN=0.5*(1.0+((2.0*RK( IM, JOIM-1)*DEEP( IM, JOIM-1))/SINH( 
* 2.*RK( IM, JOIM-1)*DEEP(IM, JOIM-1)))) 
CG( IM, JOIM-1)=C( IM, JOIM-1)*EN 
C*WILL USE SNELL’S LAW TO DETERMINE THE WAVE ANGLE HERE 
C*ANGLE OF CONTOUR WILL BE ASSUME TO BE THE SAME AS THE JMAX+1 CONTOUR 
IF(IDUM.EQ. 1)ALPH=ATAN((Y( IM, JOIM-1)-Y(IM-1,JOIM-1))/DX) 
IF (IDUM.EQ.-1)ALPH=ATAN((Y(IM+1,JOIM-1)-Y(IM, JOIM-1))/DX) 
DALPHA=ANGGEN-ALPH 
THETA( IM, JOIM-1)=ARSIN((C( IM, JOIM-1)/CGEN) *SIN(DALPHA) ) 
THETA( IM, JOIM-1)=THETA( IM, JOIM-1)+ALPH 
4 JSIM=JMAX- 4 
6 AA=0.5*(Y(IM, JSIM)+(Y(IM, USIM+1))) 
IF(AA.LT.YBAR) GO TO 8 
JSIM=JSIM-1 
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123600 


123700 
123800 


123900 
124000 
124100 
124200 
124300 
124400 
124500 
124600 
124700 
124800 
124900 
125000 
125100 
125200 
125300 
125400 
125500 
125600 
125700 
125800 
125900 
126000 
126100 
126200 
126300 
126400 
126500 
126600 
126700 
126800 
126900 
127000 
127100 
127200 
127300 
127400 
127500 
127600 
127700 
127800 
127900 
128000 
428100 
128200 
128300 
128400 
128500 
128600 
128700 
128800 
128900 
129000 
129100 
129200 
129300 
129400 
129500 
129600 
129700 
129800 
129900 
130000 
130100 
130200 
130300 
130400 
130500 
130600 
130700 
130800 


C*IF JSIM=0,THERE IS NO ADJ PT, SUB REFRAC CAN HANDLE IT. 
IF(JUSIM.EQ.O) GO TO 8 
GO TO 6 


8 RETURN 
END 


CK KK Fe oe ee kee te keke oo ke oe kek ke kk ok oko kk oo ko kok aK ok I a ok oo aK a aK ok ok ok ok ok 


SUBROUTINE WVNUM(DEPTH,T,RK) 
G=32.17 
EPS=0.001 
TWOPI=6.283185307 
SIGMA=TWOPI/T 
RK=TWOPI/(T*SQRT(G*DEPTH) ) 
DO 100 IT=1,20 
ARG=RK*DEPTH 
EK=(G*RK*TANH(ARG) )-(SIGMA**2) 
EKPR=G* (ARG*((SECH( ARG) )**2)+TANH(ARG) ) 
RKNEW=RK-EK/EKPR 
IF (ABS(RKNEW-RK) .LE.ABS(EPS*RKNEW) ) GO TO 120 
RK=RKNEW 
100 CONTINUE 
WRITE(6, 1000) IT,DEPTH, RK é 
1000 FORMAT(///,10X,"ITERATION FOR K FAILED TO CONVERGE AFTER" 
* ,3X,13,"ITERATION",/, "OUTPUT: DEPTH, RK",3X,2F 13.5) 
CALL EXIT 
120 RK=RKNEW 
IF(RK.GT.O.O) GO TO 140 
WRITE(6, 1020) DEPTH,RK 
1020 FORMAT(///,10X," RK IS NEG",/," OUTPUT DEPTH,RK",3X,2F13.5) 
CALL EXIT 
140 RETURN 
END 


Ce OR ee ee oe ee eo oe oo ok kk kak kk kk ok eo ok ok ke ok ko oo ok ok 


SUBROUTINE SMOOTH( THETA, IMAX,JMAX,IJET,SJETTY ,MMAX,Y) 


C*THIS WILL SMOOTH THE WAVE ANGLE FIELD TO ACCT FOR DIFF(ARTIFICIALLY) 


DIMENSION TEMP(60, 20),Y(60, 20), THETA(60, 20), IJET( 10) 
C*(MMAX+1) IS REQ’D BECAUSE M-GROINS HAVE M+1. REACHES OF SHORELINE. 
DO 10 M=1,MMAX+14 
IF(M.NE.1) GO TO 3 
ILEFT=2 
IRIGHT=IJET(1) 
GO TO 5 
3  IF(M.NE.MMAX+1) GO TO 4 
ILEFT=IJET(MMAX )+14 
IRIGHT=IMAX-1 
GO TO 5 
4 ILEFT=IJET(M-1)+14 
IRIGHT=IJET(M) 
5 CONTINUE 
DO 1 J=1,JMAX-1 
DO 1 I=ILEFT,IRIGHT 
IF(I.NE.ILEFT.AND.I.NE.IRIGHT) GO TO 15 
C*TO GET HERE, MUST BE ON BOUN OR ADJ TO A STRUCTURE. 
IF(I.EQ.2.0R.1.EQ. IMAX-1) GO TO 15 
C*TO GET HERE,ADJ TO A STRUCT AND CAN BE ILEFT OR IRIGHT. 
IF(Y(1,JU).GE.SJETTY) GO TO 15 
C*IF HERE, WITHIN JETTY AND ADJ TO EITHER SIDE. 
IF(I.EQ.ILEFT)TEMP(1I,JU)=0.5*(THETA(I,JU)+THETA(I+1,U)) 
IF(I.EQ.IRIGHT)TEMP(I,JU)=0.5*(THETA(I,JU)+THETA(I-1,U)) 
GO TO 1 
145 TEMP(I,J)=0.25*THETA(I-1,JU)+O.50*THETA(I,J)+O.25*THETA(I+1,U) 
4 CONTINUE 
10 CONTINUE 
DO 2 J=1,JUMAX-1 
DO 2 I=2,IMAX-1 
2 THETA(I,J)=TEMP(I,uJ) 
RETURN 
END 


CE OI IOI FO IO IO I III III I II I EI IC Oa ae ak ae 


FUNCTION SECH(A) 
SECH=1.0/COSH(A) 
RETURN 
END 
C****HERE IS WHERE THE IMSL ROUTINES MUST GO! 
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APPENDIX C 
CONTOURS AND SCHEMATIC ILLUSTRATIONS 
This appendix presents tables of the original contours at Oregon Inlet and 
the final contours for the eight numerical simulations (Tables C-1 to C-9). 


Also included are schematic illustrations of sediment volumes transported 
from the nourished region (Figs. C-1 to C-8). 
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17-ft level 


240,352 ya? 


237,284 ya? 
a —— 


3 3 
530,160 yd 522,784 yd 
SP 
4-ft level _ _ 
Ke — — Ke eee a aaa ae eee eee 7T 
Case 2.a. 
Period Considered: Twelve months, January through December, using 1975 
WIS wave hindcasts 
Sediment Budget Summary : 
Amount of sediment added: None 
c Op aA gy 
Amount of sediment transported shoreward from nourished region: 992 yd 
Amount of sediment transported seaward from nourished region: 96 vale 
3 


Net amount of sediment transported alongshore from nourished regioml9,444 yd 


Total amount of sediment transported from nourished region: 29,356 yd 


Figure C-1. Schematic illustration of sediment volumes transported 
from region, case 2.a. 
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14-ft level 
4,936 wae 
— 
1ll-ft level 
12,020 ya? 
Se 
Ke — - —- —~ —- — — -— b ——— 7-ft level —--4J 
E 1,624 ya? - 
ee aes tae ak acc FE nerf ai ee Snel SiN Lae eee | ka ee 4-ft level ea 


Case No. 2b. 
Period considered: Twelve months, January through December, using 1975 


WIS wave hindcasts, but wave angle always set equal 
to. O°. 


Sediment Budget Summary : 


Amount of sediment added: None 
Amount of sediment transported shoreward from nourished region: 1,624 ae 
Amount of sediment transported seaward from nourished region: 132 ve 

3} 


Net amount of sediment transported alongshore from nourished region-15,904 yd 


Total amount of sediment transported from nourished region: -14,148 sail 


Figure C-2. Schematic illustration of sediment volumes transported 
from region, case 2.b. 
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Case 2.cl. 


SIS es er ot — — — —~A— ~—~SY—i«*7-ft level 
460,264 ya? 


-_—-—- — _—_-— — CO - 4-ft level 


Period considered: Twelve months, January through December, using 1975 


WIS wave hindcasts. 


Sediment Budget Summary: 


Amount of sediment added: 1,452,000 Rae (on 7- and 11-ft contours) 


Amount of sediment transported shoreward from nourished region: 


Amount of sediment transported seaward from nourished region: 


Net amount of sediment transported alongshore from nourished region: 403,944 yd 


Total amount of sediment transported from nourished region: 


Figure C-3. 


460,264 yd°(31.7pc¢h 
27,808 sale (1. 9pct) 
So7sepet! 


892,016 yd*(61.4pctl 


Schematic illustration of sediment volumes transported 
from nourished region, case 2.cl. 
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ae a ie fteleve lio 


28,920 ee 
14-ft level 


5} 
138,864 yao 336,876 yd 
—_ —- 
are eS tiotia S ae ftalevel sad 
3 
637,584 yd 
414,168 yd ¢ 
eae ——p 
SS --—-3- 
466,260 yd 
Case 2.c2. 
Period considered: Twelve months, April through March, using 1975 
WIS wave hindcasts. 
Sediment Budget Summary: 
Amount of sediment added: 1,452, 000 i (on 7- and 11-ft contours) 
Amount of sediment transported shoreward from nourished region: 466,260 vai (32.1pct) 
3 
Amount of, sediment transported seaward from nourished region: 28,920 yd (2.C pct) 


Figure C-4 


3 
Net amount of sediment transported alongshore from nourished region: 421,428 yd (29.Opct) 


3 
Total amount of sediment transported from nourished region: 916,608 yd (63.lpct) 


Schematic illustration of sediment volumes transported 
from nourished region, case 2.c2. 
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—-—  17-ft level 


23,728 vdeo 
14-ft level 


3 
147,984 yd 298,056 wae 
— ae 
= = — — = = Sse ; 
3 598,400 yd 
418,072 yd 
—$»~ —_ 
Son ui eulantatasia an siete ct 7-ft level =o" 
415,784 ya? a 
Case 2.c3. 
Period considered: Twelve months, July through June, using 1975 
WIS wave hindcasts. 
Sediment Budget Summary: 
Amount of sediment added: 1,452,000 vale (on 7= and 1l-ft contour) 
Amount of sediment transported shoreward from nourished region: 415,784 sek (28.6 pet) 
Amount of sediment transported seaward from nourished region: 23,728 vale (1.6 pet) 


5 2 3 
Net amount of sediment transported alongshore from nourished region330,400 yd™ (22.8pcr) 


Total amount of sediment transported from nourished region: 769,912 ya? (53.0pct) 


Figure C-5. Schematic illustration of sediment volumes transported 
from nourished region, case 2.c3. 
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17-ft level _. — 


SS SSS OO SI 


28,452 yd> 


14-ft level ——— 


3 
151,592 yd° TIGERS Sel 
_—_—— 
3 
613,548 yd 
SECS. SC) aoe — 
> Ia 
395,556 yd 
Case 2.c4 
Period considered: Twelve months, October through September, using 1975 
WIS wave hindcasts. 
Sediment Budget Summary: 
Amount of sediment added: 1,452,000 way (on 7- and 11-ft contours). 
3 
Amount of sediment transported shoreward form nourished region: 395,556 yd (27.2 pct) 
Amount of sediment transported seaward from nourished region: 28,452 vide (2.0 pet) 
Net amount of sediment transported alongshore from nourished region: 348,144 ya? (24.0 pet) 
Total amount of sediment transported from nourished region: M12), Loe ya? (53.2 pet) 


Figure C-6. Schematic illustration of sediment volumes transported 
from nourished region, case 2.c4. 
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4,112 ya? 


— oe ee eee ee ee — ss LJ ft) level! 
3 3 
11,240 ya LE ELT TLE TEE ANS Sho 
eo ZZ 
eer ae ee E 14-f£t level Be 3 
PE SAA 3 54,528 yd 
DID ap ALT DAT ELE CLAP LTT UTALLCTE TEE PUT TEFL TET TET 
Ss ee —_— =— =— -— LI=f£t level 
32,164 ya? 
Case 3. 
Period considered: Four months, January through April, using 1975 
WIS wave hindcasts. 
Sediment Budget Summary: 
Amount of sediment added 363,000 de (on ll- and 14-ft contours). 
Amount of sediment transported shoreward from nourished region: 32,164 vde (8.9pct) 
Amount of sediment transported seaward from nourished region: 4,112 ee (1.1 pet) 


Net amount of sediment transported alongshore from nourished region: 43,708 vada (12.0 pct) 


Total amount of sediment transported from nourished region: 79,984 ya> (22.0 pct) 


Figure C-7. Schematic illustration of sediment volumes transported 
from nourished region, case 3. 
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Case 4. 275,796 yd? 


Period considered: Twelve months, January through December, using 1975 


WIS wave hindcasts. 


Sediment Budget Summary: 


Figure C-8. Schematic illustration of sediment volumes transported 


Amount of sediment added: 1,452,000 sale (on 7-, 8-, 9-, and 10-ft contours). 
Amount of sediment transported shoreward from nourished region: 275,796 ya? (19.Opct) 
Amount of sediment transported seaward from nourished region: 392,000 yaa (27.C pet) 


Net amount of sediment transported alongshore from nourished region: 96,679 vain ( 6.7pcet) 


: 3 3 
Total amount of sediment transported from nourished region: 764,475 Y¢ (52.6 pct) 


from nourished region, case 4. 
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APPENDIX D 


METHODOLOGY AND PROGRAM LISTING OF COMPUTER PROGRAM WHICH 
CONVERTS BATHYMETRIC DATA INTO MONOTONICALLY DECREASING DEPTH CONTOURS 


In order to simulate prototype shorelines (and in this case to help 
verify the numerical model via Channel Islands Harbor data), the (x, y, z) 
data points must be transformed into a form suitable for use in the model 
(i.e., bars can not be present). First, the bathymetric data have to be put 
into a form with fixed longshore and offshore spacings (i.e., ax and ay equal 
constants). This can be accomplished using one of the many available canned 
programs which do the interpolation. The problem is then one of finding the 
most suitable value of the constant, A, in the equation h = Ay2/3, 

However, as iS usually the case, the exact location of the shoreline (h = 0) 
is unknown. In addition, one requires the added constraint is required that 
the volumes of sediment (or conversely, the water above the profiles) 
balance. The problem is solved using LaGrange Multipliers and the Newton 
Raphson technique for non linear equations. 


The equation to be minimized is 


IMAX IMAX 2 


F(A,ydel,, ydel,, ... ydel Ses anh - jh ) 
3 1°? ’ (D-1) 
2 IMAX fail) ail meas; 3 pred; ; 


where A is the scale parameter in the equilibrium beach profile, ydel. are the 
locations of the shoreline for the IMAX profiles, h is the interpdlated 
depth from the survey, and Nored is the depth predicted by the equation 


2/3 


- ydel,) (D-2) 


g(A,ydel,, AAA ydel ray) = V oe Brine of i Faty - yaet ,)2/° dy 
P i=1 ydel, 
IMAX 
a 3 5/ 3am 
= 2 R Ax Alye - ydel, ) = Vane (D-3) 


where Vpreq is the predicted volume of water above the profile to the 
reference datum, Vmeas is the measured volume computed from the survey, Ax 
is the longshore distance between onshore-offshore profiles, and yf is the 
distance offshore to the last point on each of the measured profiles (it was 
a constant after the interpolation routine was used). 


102 


LaGrange Multipliers procedure says to form the quantify F* as 


eg eye (D-4) 
take the total differential of equation (D-4) 


S _ (dF dF 
dF* = dF - x dg = € dA + dtydet, J d(ydel , ) aa O56 dtydeTs aay aC) 


d d d 
ce) (4 dA + atyder 7 d(ydel, ) ieiaetie\(s cee sy) 


Rearranging 


Se eee dF d 
Ole 5 ve i) a ea mA tye, AT Se oe 


(D-6) 


It is clear that the terms in brackets in equation (D-6) must individually 
equal zero, however this leaves (IMAX + 2) unknown (udel i = to IMAX, A, and 
a4) and only (IMAX = 1) Equations. The (IMAX + 2)th equation jis taken as 
equation (D-3). The following system of equation then results: 


IMAX JMAX 
ede dam 2/3 is 2/3 
CRaecivat dat fr” oe eee ec con ee Nye ly del ys My ven er 
Vel jel Vs 
IMAX 
SNE zB Ax (Y~ - ydel rae (D-7-1) 
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wen 
Ona aydel,) ren atyde1,)= 2( Peas, j . AY > 


Bb) 


2/3) 


-1/3 2/3 


*(2/3 Aly, ar ydel,) +r ax A (ye - ydel,) 


(D-7-2) 


JMAX 278 
dF 

= » Fost - Aly .-ydel ) 

d{ydeliyay) — dlydeT yay =1 [zen meastuay. j IMAX, j IMAX 


“1/35 Als 


* (2/3 A(Y qx 57-¥92) ray) + ax A (Ye - ydel ray) 


(D-7-(IMA¥+1)) 


IMAX 


V See 


5/3 
meas ~ , (3/5 ax Aly e- ydel,) ) 


1 (D-7-(IMAX+2) ) 


Because Equations (D-7) is a system of nonlinear equations, it can not 
be written in matrix form as a [D] [x] = [E] system of equations (the 
brackets denote matrices). To solve the equations, a Newton-Raphson Iteration 
technique for nonlinear equations was used. This is done by differentiating 
each of the (IMAX + 2) equations with respect to each of the unknowns, the 
resulting equations are then linear in terms of Aa, Aydelj, . . 

Aydeltmax, 4X - The resulting matrix is inverted to obtain the a(unknown) 
and the quantities are added to the original estimates to produce a better 
estimate. This iterative procedure is continued until the changes become 
acceptably small. The solution converged rapidly. Benenaillye the first row 
of the matrix to be inverted is (aj] represents the kth row and the 1th 

column of the matrix). 


IMAX JMAX 
aq) % cS 2(y. ydel ye 
i=l jel 
JMAX 
. 4 -1/3 f ‘ 2/3 
oa? as 3 (Y. f ydel,) (reas ij any 5 ydel) ) 
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a =a (y . - ydel (h 
1, IMAX+1 3 ‘7 IMAX, j IMAX meaSimay  j 


jel 
2A(Y ray TLD crepes) 
IMAX 
ay IMAX+2 =. 05 Ls AX (Ye - yde1 )°/? ] (D-8) 
j=l 
The second row of the matrix is as follows: 
JMAX 
é 4 Ss) 6 1/3 
ao 4 Si [3 h meee i - ydel,) ie Aly, jydel,) / ] 
is VJ D 
+ AX (Ye - yde1,)°/° 
JMAX 
i 4 SAS a2 =21/3 
- » (2/3) ax A ly, ~ yder,) 71/9 
e23t= 0 
49 IMAX+1 = 2 
es 2/3 i 
a> IMAX+2 7 ox A (ye- ydel,) (D-9) 


The third row is simply these elements repeated except that the ones on the 
right-hand side of the first and last elements are changed to twos, and the 
a3 3 element is similar to the ap 9 except the ones on the right hand 

side become twos. The remaining Column elements (i.e., those when the k = 1) 
are zeroes. This process is continued to fill the array, except for the last 


row. 


The (IMAX+2)th row is as follows: 
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4IMAX+2,1 7 | Bax (Ye 


AIMAX+2,2 = -ax A (y~ - ydel,) 
AIMax+2, IMAX+1 ~ “4X A (Ve - ydelinay) 2/3 


AIMAX+2, IMAX+2 ~ © (D-i0) 


The E matrix in the [D] [x] = [LE] system of equations is 


IMAX JMAX 
28 2/3 
Se) We Bieaeth = A(yi> 2 -ydell. im wily. aasy.del) 
1 il jel meas; 4 lied 1 Is 1 
IMAX 
SUG! OF (3) Ax (Ye - ydel )°/7] 
j=l 
JMAX 
: 2/3 ane -1/3) 
ON ha ee 2(h meas, Aly, 4° ydel,) )((3) A Whe ydel, ) 
ti enAxaA (Ye - yde1,)2/)) 
JMAX 2/3 
E = - ff) 5 2(h -Aly . - ydel ) ) 
IMAX+1 aa meas Tmay  j IMAX,j IMAX 
*((§) A (yy - ydel,) nie) + (ax A (Ye - ydel,)°/°)] 
IMAX 
3 5/3 
E =e pe ese Cm) MAA Gye mvdelll =) mand) i= av ] 
IMAX+2 72 5 f I meas (p-11) 
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The [D] [x] = [E] system of equations was then solved, as explained 
previously, by solving the x column vector (which represents the changes in 
the unknowns, AA, Aydel] ... Aydelymax. AA), adding these changes to the 
respective variables and iterating until a final solution is obtained. 


The computer program which did these calculations for the Channel Island 


Harbor simulation follows. A user-supplied matrix inversion routine is 
required (Line 37,200). 
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GRESET FREE 
CreceeercoveessPROGRAM CIH/BVALUE 1 
FILE S(KIND=PACK, TITLE="CIH42076A" ,FILETYPE=7) 
FILE 6(KIND=REMOTE ) 
C*THIS PROGRAM USES THE INTERPOLATED PROFILES OF CIH. 
C*IT FINDS THE LOCATION OF THE SHORELINE, YDEL AND THE BEST 
C*FIT LEAST SQUARES "“B" VALUE FOR H=BY**2/3 
C*USES LAGRANGE MULTIPLIERS TO CONSTRAIN THE VOLUMES(SO THEY ARE EQUAL) 
C*THEN IT USES NEWTON-RAPHSON ITER FOR NON-LIN EQS 
DIMENSION X(40) 
DIMENSION WKAREA(600) , AMATRX(23,23),BMATRX(23, 1) 
DIMENSION Y(40,20),Z(40, 20), YDEL(40), JBEGIN( 40), YDELI (40) 
DIMENSION DYTWO( 40, 20) ,DYONE (40, 20) ,.DYMTWO( 40, 20) ,.DYMONE (40, 20) 
NIMENSION DYMFOR( 40,20) ,0YFOR( 40,20). YDONE( 40,20), YOMTWO( 40, 20) 
DIMENSION YOMONE (40,20), YETWO( 40) , YEONE (40) , YEMONE (40) 
DIMENSION YEMTWO( 40), YEMFOR( 40), YEFIVE( 40) 
EXPON=2./3. 
THIRD=0. 3333333333333333 
C*FIRST READ IN THE PROFILES FROM DISKPACK. 
DO 1 I=1,34 
DO 1 J=1,15 
1 READ(5, 100) X(I).Y¥(1,U),2Z(1,u) 
100 FORMAT(14X,F6.0,F5.0,F5.0) 
C*NOW WE MUST GET A FIRST APPROX FOR YDEL 
C*WE WILL USE LINEAR INTERPOLATION TO DETERMINE IT. 
IBEGIN=1 
IMAX=21 
JUMAX=15 
C*CHANGE PROFILE TO SPAN 1 TO IMAX(IF ALREADY DONE,WON’T HARM THINGS) 
ITEMP1=1 
ITEMP2=IMAX-IBEGIN+1 
K=-1 
OO 777 I=1,I1TEMP2 
K=K+1 
DO 777 J=1,JUMAX 
Y(1I,U)#Y(IBEGIN+K,J) 
777 2(1,.JU)=Z(IBEGIN+K, J) 
IMAX=I TEMP2 
OX= 100.00 
DO 2 I=1, IMAX 
DO 3 J=1,JUMAX 
IF(Z(1,JU).GE.0.0) GO TO 3 
C*FIRST NEG POINT ON THE PROFILE IS SEAWARD OF 220.0 
C* WE MUST ALSO REMEMBER THIS LOCATION. 
C*IF Z(1,1)<O.,CHOOSE ARBITRARY PT, ROUTINE ITERATES TO SOLN. 
ZDUM=1.0 
IF(JU.NE. 1) ZDUM2Z(1,JU-1) 
YDUM=Y(I,JU)-50.0 
IF(U.NE.1) YOUM2Y(I,uU-1) 
DELY=ZDUM/((ZDUM-Z(1,JU))/(Y¥(1,JU)-YDUM) ) 
YDEL(I )=YDUM+DELY 
JBEGIN(I)=J 
GO TO 2 
3 CONTINUE 
2 CONTINUE 
C*THE VALUES FOR Z ARE NEG ON FILE, MUST NOW MAKE POS. 
C*THE Z VALUES ARE ALSO *10. 
OO 35 I=1, IMAX 
OO 35 J=JBEGIN(1).UMAX 
35 2(1,JU)=-Z(1,u)/10.0 
C*MUST INITIALIZE “B" SO WILL MAKE A FIRST GUESS. 
C*MUST ALSO GUESS LAMBDA (XLAMB) 
B=0.30 
XLAMB=-2.0 
00 10 ITER=1,100 
C*LET’S CALCULATE THE VOL OF WATER ABOVE THE PROFILE,VMEAS. 
C*ITS OUR CONSTRAINT,BUT SINCE YDEL IS NOT KNOWN,A PRIORI,IT WILL CHANGE 
VMEAS=0.0 
DO 200 1=1, IMAX 
00 200 J=JBEGIN(I), UMAX 
I1F(JU.NE. JBEGIN(I)) GO TO 201 
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VME SS=VMEAS*OX*2(1.U)9(0.5*(V(1,u)+¢(1,U+1))-YDEL(1)) 
GO TO 200 
201 IF(J.EQ. UMAX) GO TO 202 
VMEAS=VMEAS+DX*0.5*(V¥(1,u+1)-Y¥(I.U-1))*Z2(1.uU) 
GO TO 200 
202 VMEAS=VMEAS+DX*2(1,U)*(Y¥(I,u)-0.5*(¥(J.U)+¥(1,U-1))) 
200 CONTINUE 
C*PRIOR TO EQS,COMPUTE AND STORE SEVERAL VALUES WE NEED OVER AND OVER 
C*BECAUSE COMPUTER CAN’T RAISE A NEG VALUE TO AN EXPONENT 
C*MUST PRESERVE THE SIGN. 
DO 400 II=1, IMAX 
DO 401 JJ=JBEGIN(II),UMAX 
ARG1=¥(11,JJ)-YDEL(II) 
DYSIGN=SIGN(1.,ARG1) 
OY=ABS(Y(II,JUJ)-YDEL(II)) 
DYTWO(I1,JJU)=DY**EXPON 
OVONE (II ,JJU)=DYSIGN*DY** THIRD 
OVMTWO( II, JJ) =DY**(-EXPON) 
OYMONE (II , JJ) =DYSIGN*DY * *(-THIRD) 
DYMFOR(II ,JJU)=DY**(-2.*EXPON) 
OYFOR(II, JJ) =DY**(2. *EXPON) 
401 CONTINUE 
ARG2=1400. -YDEL(II) 
DSIGN=SIGN(1. ,ARG2) 
DYE=ABS(ARG2) 
YETWO( II )=O0YE**EXPON 
YEONE(I1I)=DSIGN*DYE**THIRD 
YEMONE (II )=DOSIGN*DYE**(-THIRD) 
YEMTWO( II) =DYE**(-EXPON) 
YEMFOR(II)=DYE**(-2. *EXPON) 
YEFIVE(II)=DSIGN*DYE**(5.*THIRD) 
400 CONTINUE 
C*LET’S INPUT THE FIRST ROW OF THE MATRIX, A 
SUM1B=0.0 
DO 300 II=1,1MAX 
DO 300 JJ=UBEGIN( II), JMAX 
300 SUM1B=SUM1B+2. *DYFOR(II, Ju) 
AMATRX(1, 1)=SUM1B 
SUMLAM=0.0 
DO 305 K=1, IMAX 
SUM1K=0.0 
DO 306 JJ=JBEGIN(K) , MAX 
306 SUM1K=SUM1K+2. *EXPON*DYMONE(K, JJ) *(Z(K,UJ)-2. *B* 
* DYTWO(K,JUJ)) 
SUMLAM=SUMLAM-O. 6*DX *YEFIVE(K) 
305 AMATRX(1,K+1)=SUMIK 
AMATRX( 1, IMAX+2)=SUMLAM 
C*NOW THE MIDDLE ROWS OF THE AMATRX. 
DO 410 LROW=2,IMAX+1 
SUM2B=0.0 
IT=LROW-1 
DO 415 JJ=JUBEGIN(II),UMAX 
415 SUM2B=SUM2B+2. *EXPON*Z(II,JUJ) *DYMONE(II,JJ)-4. *EXPON* 
* B*DYONE(I1,JJ) 
AMATRX(LROW, 1) =SUM2B+XLAMB*DX*YETWO(I1) 
00 430 II=1, IMAX 
SUM2Y=0.0 
00 425 VJ=JBEGIN(II),JMAX 
425 SUM2Y=SUM2Y+2. *EXPON*THIRD*B*Z(II,JJU) *DYMFOR(II,JUJ)+THIRD*EXPON 
* *2.*B*B*DYMTWO(II, JU) 
IF((2141).EQ.LROW) GO TO 431 
AMATRX(LROW, I1I1+1)=0.0 
GO TO 430 
431 AMATRX(LROW,I1+1)=SUM2Y-XLAMB*EXPON*DX*B*YEMONE (II) 
430 CONTINUE 
410 AMATRX(LROW, IMAX+2)=DX*B*VETWO(LROW- 1) 
C*NOW THE LAST ROW OF THE MATRIX A 
SUMFB=0.0 
DO 450 Tl=1, IMAX 
450 SUMFB=SUMFB+0O.6*DX*VEFIVE(II) 
AMATRX( IMAX+2, 1) =SUMFB 
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14300 
14400 
14500 
14600 
14700 
14800 
14900 
15000 
15100 
15200 
15300 
15400 
15500 
15600 
15700 
15800 
15900 
16000 
16100 
16200 
16300 
16400 
16500 
16600 
16700 
16800 
16900 
17000 
17100 
17200 
17300 
17400 
17500 
17600 
17700 
17800 
17900 
18000 
18100 
18200 
18300 


DO 453 1I=1, IMAX 
453 AMATRX(I1MAX+2,11+1.)2-DX*B*VYETWO(II) 
AMATRX( IMAX+2, IMAX+2)=0.0 
C*NOW MUST INPUT THE BMATRX. 
SUMF 1A=0.0 
SUMF 1B=0.0 
DO 455 I1=1, IMAX 
SUMF 1B=SUMF 1B+XLAMB*O.6*OX*YEFIVE(II) 
DO 455 JJ=JBEGIN( II) , MAX 
455 SUMF 1A=SUMF 1A-2.*(Z(II,JJU)-B*DYTWO(II,JJU)) *OYTWO(II, JU) 
BMATRX(1, 1)=-(SUMF 1A-SUMF 1B) 
DO 460 I1=1, IMAX 
SUMFI1I=0.0 
DO 462 JJ=JBEGIN( II) , MAX 
462 SUMFII=SUMFII+2.*(Z(II,JUJU)-B*DYTWO(II,JJ) ) *EXPON*B*DYMONE (II, JJ) 
SUMF I 1=SUMFII+XLAMB*DX*B*YETWO( IT) 
460 BMATRX(II+1,1)=-SUMFII 
SUMV=0.0 
DO 465 IIl=1, IMAX 
465 SUMV=SUMV+0.6*DX*B*YEFIVE(II) 
BMATRX(IMAX+2, 1)=-(SUMV-VMEAS) 
C*NEXT LET’S CALL THE MATRIX INVERSION ROUTINE VIA IMSL 
CALL LEQT2F(AMATRX, 1, IMAX+2,23,BMATRX,3,WKAREA, IER) 
C*THE SOLN IS RETURNED IN THE VECTOR BMATRX 
C*FINALLY, WE MUST UPDATE THE X VECTOR IN AX=B. 
B=B+BMATRX(1,1) 
XLAMB=XLAMB+BMATRX( IMAX+2, 1) 
00 470 I1=1, IMAX 
470 YDEL(II)=YOEL(II)+BMATRX(II+1,1) 
C*CHECK THE CRITERION FOR COMPLETION 
SUMVEC=0.0 
DO 475 IIl=1, IMAX 
475 SUMVEC=SUMVEC+ABS(BMATRX(II,1)) 
IF (SUMVEC.LT. (0. 1*( IMAX+2))) GO TO 11 
WRITE(6,*/) B,ITER,(1,YOEL(I),121, IMAX),XLAMB 
10 CONTINUE 
11 CONTINUE 
C*LET’S WRITE IT ALL OUT. 
WRITE(6, */) ITER.B. (1, VOEL(I), 124, FMAX) 
STOP 
ENO 
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APPENDIX E 


USER DOCUMENTATION AND INPUT AND OUTPUT FOR PROGRAM VERIFICATION 


The computer program presented in Appendix B was run on a Burroughs 
B-7700 computer. The B7000/B6000 series FORTRAN language was designed so 
several existing programs written in FORTRAN would be compatible with minimal 
changes. It was designed to be compatible with Fortram IV, H level and to 
contain ANSI X3.9-1966 Standard FORTRAN as a subset. 


Line 37,200 of the coding (see App. B) requires a subroutine from the 
IMSL subroutine package, LEQT1B and its associated subroutines. If the 
user's computing center has access to this package of subroutine programs 
they need only bind them to the program (note: copyright laws prohibited the 
inclusion of the IMSL coding). If not, a substitute subroutine must be user 
supplied. It must facilitate the solution of a banded storage mode matrix. 


The program input will be described here using a card deck set-up, 
however, the use of diskpack or magnetic tape input follows directly. Lines 
3100, 4100, 5500, 5900, 6800, 7500, and 12,900 are read statements. The 
cards used for the simulation presented in this appendix are shown in Figure 
E-1. The first card contains the value of WDEPTH, the depth of water (in 
meters) to which the input wave conditions are to be transformed (a partial 
list of variables used in the program is presented beginning on page A-8 of 
Appendix A). The format statements are obviously in the program coding. 


The second data input card is read by line 4100 where the variables 
SJETTY, BERM, SFACE, and DIAM are required (length of the structure, berm 
height, shore face slope, and sediment diameter, respectively). 


Lines 5500 reads MMAX, the number of structures to be simulated (as 
set-up here, a maximum of 10 structures can be modeled, however, appropriate 
changes in array dimensions would allow additions (structures). Line 5900, 
which is in a "DO" loop reads the lesser I grid value adjacent to where the 
structure is desired. The number of structures, MMAX, determines the number 
of data cards required here; 3 structures require 3 cards with the 3 I grid 
locations (note, the present configuration of the refraction and diffraction 
subroutines requires evenly spaced structures, however this can be altered if 
necessary). 


The parameter ADEAN, which represents the value of A in the equilibrium 
profile used is the next value input (line 6800). As mentioned previously, 
whenever possible a site-specific value should be used. The space-step and 
time-step (DX and DELT in the coding) are input next (line 7500). 


The last input values are the wave data, HS, T, ALPWIS read by line 
12,900. This statement is in a loop made by the unconditional GO TO statement 
(line 16,400) and the read statement. There is an action specifier included 
in the read statement to transfer the program to statement 1000, thereby stop- 
ping execution of the program once all the wave climate data have been used. 
The number of data cards required for this read statement is dictated by the 
length of the simulation and the time-step used. 


The input file and output for program verification follow. 


111 


WDEPTH 
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Figure E-1. Card deck input for program verification. 
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FILE DUM 
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